{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#首先导入必要的工具包\n",
    "import numpy as np  # 矩阵操作\n",
    "import pandas as pd # SQL数据处理\n",
    "\n",
    "from sklearn.metrics import r2_score  #评价回归预测模型的性能\n",
    "\n",
    "import matplotlib.pyplot as plt   #画图\n",
    "import seaborn as sns\n",
    "\n",
    "# 图形出现在Notebook里而不是新窗口\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数据探索"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>dteday</th>\n",
       "      <th>season</th>\n",
       "      <th>yr</th>\n",
       "      <th>mnth</th>\n",
       "      <th>holiday</th>\n",
       "      <th>weekday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>weathersit</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>casual</th>\n",
       "      <th>registered</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2011-01-01</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.344167</td>\n",
       "      <td>0.363625</td>\n",
       "      <td>0.805833</td>\n",
       "      <td>0.160446</td>\n",
       "      <td>331</td>\n",
       "      <td>654</td>\n",
       "      <td>985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>2011-01-02</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.363478</td>\n",
       "      <td>0.353739</td>\n",
       "      <td>0.696087</td>\n",
       "      <td>0.248539</td>\n",
       "      <td>131</td>\n",
       "      <td>670</td>\n",
       "      <td>801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>2011-01-03</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.196364</td>\n",
       "      <td>0.189405</td>\n",
       "      <td>0.437273</td>\n",
       "      <td>0.248309</td>\n",
       "      <td>120</td>\n",
       "      <td>1229</td>\n",
       "      <td>1349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2011-01-04</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>0.212122</td>\n",
       "      <td>0.590435</td>\n",
       "      <td>0.160296</td>\n",
       "      <td>108</td>\n",
       "      <td>1454</td>\n",
       "      <td>1562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>2011-01-05</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.226957</td>\n",
       "      <td>0.229270</td>\n",
       "      <td>0.436957</td>\n",
       "      <td>0.186900</td>\n",
       "      <td>82</td>\n",
       "      <td>1518</td>\n",
       "      <td>1600</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant      dteday  season  yr  mnth  holiday  weekday  workingday  \\\n",
       "0        1  2011-01-01       1   0     1        0        6           0   \n",
       "1        2  2011-01-02       1   0     1        0        0           0   \n",
       "2        3  2011-01-03       1   0     1        0        1           1   \n",
       "3        4  2011-01-04       1   0     1        0        2           1   \n",
       "4        5  2011-01-05       1   0     1        0        3           1   \n",
       "\n",
       "   weathersit      temp     atemp       hum  windspeed  casual  registered  \\\n",
       "0           2  0.344167  0.363625  0.805833   0.160446     331         654   \n",
       "1           2  0.363478  0.353739  0.696087   0.248539     131         670   \n",
       "2           1  0.196364  0.189405  0.437273   0.248309     120        1229   \n",
       "3           1  0.200000  0.212122  0.590435   0.160296     108        1454   \n",
       "4           1  0.226957  0.229270  0.436957   0.186900      82        1518   \n",
       "\n",
       "    cnt  \n",
       "0   985  \n",
       "1   801  \n",
       "2  1349  \n",
       "3  1562  \n",
       "4  1600  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#读取数据\n",
    "# 调用csv中的数据\n",
    "data = pd.read_csv(\"day.csv\")\n",
    "\n",
    "#通过观察前5行，了解数据每列（特征）的概况\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(731, 16)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#数据基本信息（条数，特征数）\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'DataFrame' object has no attribute 'column'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-23-fe900828460d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;31m#用于后续显示权重系数对应的特征\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[0mcolumns\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mX\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcolumn\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m \u001b[0mX\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhead\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mD:\\Anaconda\\lib\\site-packages\\pandas\\core\\generic.py\u001b[0m in \u001b[0;36m__getattr__\u001b[1;34m(self, name)\u001b[0m\n\u001b[0;32m   3612\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mname\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_info_axis\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   3613\u001b[0m                 \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 3614\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mobject\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__getattribute__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   3615\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   3616\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__setattr__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: 'DataFrame' object has no attribute 'column'"
     ]
    }
   ],
   "source": [
    "# 从原始数据中分离输入特征x和输出y\n",
    "y = data['cnt'].values\n",
    "X = data.drop(['cnt','dteday','mnth','casual','registered'],axis = 1)\n",
    "\n",
    "#用于后续显示权重系数对应的特征\n",
    "columns = X.column\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "对X进行处理，除了删除cnt之外,还按照题目要求抛弃了casual和registered。除此之外，根据之前的数据探索结果，发现月份和季节高度相关，说明存在数据冗余，于是删除与cnt相关性较小的mnth"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "这边的报错通过查询之后，是.pyc文件存在问题，但是删除.pyc文件之后错误依然存在，目前无法解决该问题。之后会对DataFrame进行进一步的了解，看是否可以对问题有更深层次的认识。如果助教遇到过类似问题，希望可以给予指导，谢谢。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(365, 11)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#数据准备\n",
    "#按照作业要求将数据分割为训练数据和测试数据\n",
    "\n",
    "#将数据分割训练数据与测试数据\n",
    "\n",
    "#将2011年的数据作为训练数据\n",
    "data_train = X[data.yr == 0]\n",
    "\n",
    "#将2012年的数据作为测试数据\n",
    "\n",
    "data_test = X[data.yr == 1]\n",
    "#查看训练数据条数\n",
    "data_train.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(366, 11)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#查看测试数据条数\n",
    "data_test.shape "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(366,)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#对X_train, X_test, y_train, y_test进行赋值\n",
    "X_train = data_train\n",
    "X_test = data_test\n",
    "\n",
    "#对cnt进行测试数据和训练数据的划分\n",
    "y_train = pd.Series(y[0:365])\n",
    "y_test = pd.Series(y[365:len(y)])\n",
    "#X_train.head()\n",
    "#X_test.head()\n",
    "y_test.shape"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "之前在划分数据的时候曾使用data_train = data[data.yr == 0]来划分，但是并不会与X进行整合，所以只能采用相对简单的办法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda\\lib\\site-packages\\ipykernel_launcher.py:20: FutureWarning: reshape is deprecated and will raise in a subsequent release. Please use .values.reshape(...) instead\n",
      "D:\\Anaconda\\lib\\site-packages\\sklearn\\utils\\validation.py:475: DataConversionWarning: Data with input dtype int64 was converted to float64 by StandardScaler.\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "D:\\Anaconda\\lib\\site-packages\\ipykernel_launcher.py:21: FutureWarning: reshape is deprecated and will raise in a subsequent release. Please use .values.reshape(...) instead\n"
     ]
    }
   ],
   "source": [
    "#由之前的数据分析可以看出，某些特征存在特征差异较大的情况，所以进行数据标准化预处理\n",
    "#标准化的目的在于避免原始特征值差异过大，导致训练得到的参数权重不归一，无法比较各特征的重要性\n",
    "\n",
    "#X_train = data_train\n",
    "#X_test = data_test\n",
    "\n",
    "# 数据标准化\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# 分别初始化对特征和目标值的标准化器\n",
    "ss_X = StandardScaler()\n",
    "ss_y = StandardScaler()\n",
    "\n",
    "# 分别对训练和测试数据的特征以及目标值进行标准化处理\n",
    "X_train = ss_X.fit_transform(X_train)\n",
    "X_test = ss_X.transform(X_test)\n",
    "\n",
    "#对y做标准化不是必须\n",
    "#对y标准化的好处是不同问题的w差异不太大，同时正则参数的范围也有限\n",
    "y_train = ss_y.fit_transform(y_train.reshape(-1, 1))\n",
    "y_test = ss_y.transform(y_test.reshape(-1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'fs' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-21-d4451c3ac2ad>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     16\u001b[0m \u001b[1;31m# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     17\u001b[0m \u001b[1;31m#fs = pd.DataFrame({\"columns\":list(columns), \"coef\":list((lr.coef_.T))})\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 18\u001b[1;33m \u001b[0mfs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msort_values\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mby\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'coef'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mascending\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'fs' is not defined"
     ]
    }
   ],
   "source": [
    "# 线性回归\n",
    "#class sklearn.linear_model.LinearRegression(fit_intercept（拟合截距项）=True, normalize=False（因为已经对数据进行了标准化）, copy_X=True, n_jobs=1)\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# 使用默认配置初始化\n",
    "lr = LinearRegression()\n",
    "\n",
    "# 训练模型参数\n",
    "lr.fit(X_train, y_train)\n",
    "\n",
    "# 预测\n",
    "y_test_pred_lr = lr.predict(X_test)\n",
    "y_train_pred_lr = lr.predict(X_train)\n",
    "\n",
    "\n",
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef\":list((lr.coef_.T))})\n",
    "fs.sort_values(by=['coef'],ascending=False)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "因为前面的columns赋值存在问题，所以该处不能执行打印权重系数操作,因此不能从权重系数中查看对应特征的重要性。（该处绝对值越大说明越重要）"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "下面为模型评价"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of LinearRegression on test is -0.7551751016700898\n",
      "The r2 score of LinearRegression on train is 0.7585428339129555\n"
     ]
    }
   ],
   "source": [
    "# 使用r2_score评价模型在测试集和训练集上的性能，并输出评估结果\n",
    "#测试集\n",
    "print ('The r2 score of LinearRegression on test is', r2_score(y_test, y_test_pred_lr))\n",
    "#训练集\n",
    "print ('The r2 score of LinearRegression on train is', r2_score(y_train, y_train_pred_lr))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "该处出现负值不知道错误出现在哪里，通过查询文档，出现负值被认为是可能出现的情况（在模型选择比较差的情况下）。\n",
    "Best possible score is 1.0 and it can be negative (because the model can be arbitrarily worse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFsCAYAAAADhPr/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHjBJREFUeJzt3XuUXGWZ7/HvQ9IShkAw0IFAjAEE5Z5gExMRjMYLI0hwjYoc5LK4REFm4XUUPQ4to6McUJcc0TPxBgjKbUQYlTNwIpkIR2ESJtwMmqhwTIhJiHIJCiThOX/UTmw63elKd1X3m67vZ61aXbX3u/d+3tqd/Orde9fuyEwkSVKZthvqAiRJUu8MakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtVpeRDwUETOGuo6hFBHviIjfR8TaiJgyiNtdGxH79DLv9Ii4s0HbeSQi3tSIdUmDzaDWsNbTf9DdAyAzD8rMeX2sZ1JEZESMbFKpQ+1S4LzMHJ2Z/9V9ZtX3Z6pgXR4RX4qIEQPdaLW93w50PdJwZlBLBSjgA8DLgYf6aHNYZo4GXg+cCJzR9KokGdRS11F3REyNiAUR8VRErIyIL1XN5lc/n6hGldMjYruI+O8R8WhErIqIqyJiTJf1nlrNWxMRn+62nc6IuDEiro6Ip4DTq23/PCKeiIgVEfHViHhJl/VlRJwbEUsi4umI+KeI2Lda5qmIuL5r+2597LHWiNg+ItYCI4D7IuI3fb1fmbkUuAuY3GX9YyLiW1XdyyPisxtH3BHxioj4j4h4MiIej4jruvXpFdXzXSPilqov9wD7dmm32RGNiJgXEWdVz/eNiJ9W7/XjEXFNROzSy3vR2z6WimRQSy/2FeArmbkztaC4vpp+dPVzl+pw7c+B06vHG4B9gNHAVwEi4kDga8DJwHhgDLBXt23NAm4EdgGuATYAHwJ2A6YDM4Fzuy1zDPBqYBrwD8CcahsvAw4GTuqlXz3WmpnPVaNkqI2Y9+158b+KiFcBRwFLu0y+ElgPvAKYArwFOKua90/AbcBLgQnA/+xl1ZcDz1J7v85g60bsAXwe2BM4gNr70dlL2972sVQkg1qt4IfVKPWJiHiCWoD2Zh3wiojYLTPXZuYvttD2ZOBLmfnbzFwLXAC8pxr1vRP4t8y8MzOfB/4R6H5j/Z9n5g8z84XM/EtmLszMX2Tm+sx8BPgXaoeZu7o4M5/KzIeAB4Hbqu0/CdxKLSS3ttZ63RsRzwCLgXlU72NE7A78LfDBzHwmM1cBXwbeUy23jtqh9T0z89nM3OwCsWr0/XfAP1breJBa+NclM5dm5u3VB4/VwJfY/L3baGv2sTTkDGq1ghMyc5eNDzYfpXZ1JrA/8HBE/GdEHLeFtnsCj3Z5/SgwEti9mvf7jTMy88/Amm7L/77ri4jYPyJ+FBF/qA6H/zO10XVXK7s8/0sPr0fTsy3VWq/Dq/WfCLwG2LGa/nKgDVjR5cPQvwDjqvn/QG3Ee0/UrrDvaaTcXtXT9T15tId2PYqIcRFxbXXY/SngajZ/7zbamn0sDTmDWuoiM5dk5knUQuZi4MaI2JHNR8MAj1ELqY0mUjv8uxJYQe0wLwARsQOwa/fNdXv9deBhYL/qsOwnqQVcI2yp1rplzfXAz6kdJYBauD4H7NblA9HOmXlQtcwfMvPszNwTeB/wtY3npbtYXdXzsm41bvRM9fNvukzbo8vzz1N7Pw+t3rv30st7t4V9LBXJoJa6iIj3RkR7Zr4APFFN3kAtSF6gdn53o+8DH4qIvSNiNLUR8HWZuZ7auee3R8Rrqwu8PkPfobsT8BSwtjoPfE7DOrblWvvjC8DsiNgjM1dQOwf9xYjYubpwbd+IeD1ARLwrIjZ+aPkTtUDd0HVlmbkB+AHQGRF/U53jP63L/NXAcuC9ETGiGpV3PZ++E7CW2sV+ewEf663wLexjqUgGtfRixwAPVVdCfwV4T3Ve9c/A54C7qsO704BvA9+ldkX476hdCPX3ANU55L8HrqU2un4aWEVt5NmbjwL/rWr7DeC6LbTdWr3W2h+Z+QDwH/w1EE8FXgL8kloY30jtojCAI4C7q/f0FuD8zPxdD6s9j9qh9T8AVwDf6Tb/7Gp7a4CDgP/bZd5nqB2afxL4MbXQ702P+3jLPZaGTmT2dERPUiNVo9gnqB3W7imkJKlHjqilJomIt1eHcXekduevB4BHhrYqSdsag1pqnlnULuJ6DNiP2iFWD2FJ2ioe+pYkqWCOqCVJKtig/iGA3XbbLSdNmjSYm5QkqTgLFy58PDPb62k7qEE9adIkFixYMJiblCSpOBFR9533PPQtSVLBDGpJkgpmUEuSVLBBPUctSerZunXrWLZsGc8+691Mh5NRo0YxYcIE2tra+r0Og1qSCrBs2TJ22mknJk2aRESj/miahlJmsmbNGpYtW8bee+/d7/V46FuSCvDss8+y6667GtLDSESw6667DvgoiUEtSYUwpIefRuxTg1qSpIJ5jlqSCtTZOfjrGzFiBIcccgjr169n77335rvf/S677LLLVm/rrLPO4sMf/jAHHnjgi6ZfccUVLFiwgK9+9atbvU6A0aNHs3bt2rrazpgxg0svvZSOjo5N0xYsWMBVV13FZZdd1q/tDxVH1JIkAHbYYQcWLVrEgw8+yNixY7n88sv7tZ5vfvObm4V0CTo6Opoe0hs2bGj4OvsM6ogYFRH3RMR9EfFQRHymmn5FRPwuIhZVj8kNr06SNCSmT5/O8uXLN72+5JJLOOKIIzj00EO58MILAXjmmWc49thjOeywwzj44IO57rrrgNpoduPtor/zne+w//778/rXv5677rpr0/pOP/10brzxxk2vR48eDcDatWuZOXMmhx9+OIcccgg333zzZrWtWLGCo48+msmTJ3PwwQfzs5/9rK4+zZs3j+OOOw6Azs5OzjjjDGbMmME+++zzogC/+uqrmTp1KpMnT+Z973vfpvA955xz6Ojo4KCDDtr0HkDt9tgXXXQRr3vd67jhhhvqqmVr1HPo+zngjZm5NiLagDsj4tZq3scy88YtLCtJ2sZs2LCBuXPncuaZZwJw2223sWTJEu655x4yk+OPP5758+ezevVq9txzT3784x8D8OSTT75oPStWrODCCy9k4cKFjBkzhje84Q1MmTJli9seNWoUN910EzvvvDOPP/4406ZN4/jjj3/RRVnf+973eOtb38qnPvUpNmzYwJ///Od+9fPhhx/mjjvu4Omnn+aVr3wl55xzDkuXLuW6667jrrvuoq2tjXPPPZdrrrmGU089lc997nOMHTuWDRs2MHPmTO6//34OPfTQTXXfeeed/aqjL30GdfWH7jeeFGirHv4Ra0kaZv7yl78wefJkHnnkEV796lfz5je/GagF9W233bYpZNeuXcuSJUs46qij+OhHP8rHP/5xjjvuOI466qgXre/uu+9mxowZtLfX/kjUiSeeyK9//est1pCZfPKTn2T+/Plst912LF++nJUrV7LHHntsanPEEUdwxhlnsG7dOk444QQmT+7fAd1jjz2W7bffnu23355x48axcuVK5s6dy8KFCzniiCM2vSfjxo0D4Prrr2fOnDmsX7+eFStW8Mtf/nJTUJ944on9qqEedZ2jjogREbEIWAXcnpl3V7M+FxH3R8SXI2L7plUpSWq6jeeoH330UZ5//vlN56gzkwsuuIBFixaxaNEili5dyplnnsn+++/PwoULOeSQQ7jgggu46KKLNltnb19PGjlyJC+88MKm9T///PMAXHPNNaxevZqFCxeyaNEidt99982+h3z00Uczf/589tprL0455RSuuuqqfvV3++3/GlsjRoxg/fr1ZCannXbapr7+6le/orOzk9/97ndceumlzJ07l/vvv59jjz32RXXtuOOO/aqhHnVd9Z2ZG4DJEbELcFNEHAxcAPwBeAkwB/g4sNleiojZwGyAiRMnNqhsSduieq9kbvQVz9o6Y8aM4bLLLmPWrFmcc845vPWtb+XTn/40J598MqNHj2b58uW0tbWxfv16xo4dy3vf+15Gjx7NFVdc8aL1vOY1r+H8889nzZo17Lzzztxwww0cdthhQO287sKFC3n3u9/NzTffzLp164Da4fNx48bR1tbGHXfcwaOPbv7XIB999FH22msvzj77bJ555hnuvfdeTj311Ib0febMmcyaNYsPfehDjBs3jj/+8Y88/fTTPPXUU+y4446MGTOGlStXcuuttzJjxoyGbLMvW/X1rMx8IiLmAcdk5qXV5Oci4jvAR3tZZg61IKejo8ND5pJUh6H+sDJlyhQOO+wwrr32Wk455RQWL17M9OnTgdqFX1dffTVLly7lYx/7GNtttx1tbW18/etff9E6xo8fT2dnJ9OnT2f8+PEcfvjhmy7MOvvss5k1axZTp05l5syZm0akJ598Mm9/+9vp6Ohg8uTJvOpVr9qstnnz5nHJJZfQ1tbG6NGjex1RH3vssZvusT19+nQ+8IEP9NnvAw88kM9+9rO85S1v4YUXXqCtrY3LL7+cadOmMWXKFA466CD22WcfjjzyyPrfzAGK2inoLTSIaAfWVSG9A3AbcDGwMDNXRO24xpeBZzPzE1taV0dHR268ElBS63FE3bvFixdzwAEHDHUZaoKe9m1ELMzMjl4WeZF6RtTjgSsjYgS1c9rXZ+aPIuKnVYgHsAh4/9aVLkmS+lLPVd/3A5tdT5+Zb2xKRZIkaRPvTCZJhejrVKS2PY3Ypwa1JBVg1KhRrFmzxrAeRjb+PepRo0YNaD3+UQ5JKsCECRNYtmwZq1evHupS1ECjRo1iwoQJA1qHQS1JBWhra2Pvvfce6jJUIA99S5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBTOoJUkqmEEtSVLBDGpJkgpmUEuSVDCDWpKkghnUkiQVzKCWJKlgBrUkSQUzqCVJKphBLUlSwQxqSZIKZlBLklQwg1qSpIIZ1JIkFcygliSpYAa1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBTOoJUkqWJ9BHRGjIuKeiLgvIh6KiM9U0/eOiLsjYklEXBcRL2l+uZIktZZ6RtTPAW/MzMOAycAxETENuBj4cmbuB/wJOLN5ZUqS1Jr6DOqsWVu9bKseCbwRuLGafiVwQlMqlCSphdV1jjoiRkTEImAVcDvwG+CJzFxfNVkG7NWcEiVJal11BXVmbsjMycAEYCpwQE/Nelo2ImZHxIKIWLB69er+VypJUgvaqqu+M/MJYB4wDdglIkZWsyYAj/WyzJzM7MjMjvb29oHUKklSy6nnqu/2iNiler4D8CZgMXAH8M6q2WnAzc0qUpKkVjWy7yaMB66MiBHUgv36zPxRRPwSuDYiPgv8F/CtJtYpSVJL6jOoM/N+YEoP039L7Xy1JElqEu9MJklSwQxqSZIKZlBLklSwei4mk6QidXY2tp1UIkfUkiQVzKCWJKlgBrUkSQUzqCVJKphBLUlSwQxqSZIKZlBLklQwg1qSpIIZ1JIkFcygliSpYAa1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBVs5FAXIGlwdXY2tl0zDOW2pdI4opYkqWAGtSRJBTOoJUkqmEEtSVLBDGpJkgpmUEuSVDCDWpKkghnUkiQVzKCWJKlgBrUkSQUzqCVJKphBLUlSwQxqSZIK1mdQR8TLIuKOiFgcEQ9FxPnV9M6IWB4Ri6rH25pfriRJraWeP3O5HvhIZt4bETsBCyPi9mrelzPz0uaVJ0lSa+szqDNzBbCiev50RCwG9mp2YZIkqb4R9SYRMQmYAtwNHAmcFxGnAguojbr/1MMys4HZABMnThxguZK09To7m9NWGgx1X0wWEaOBfwU+mJlPAV8H9gUmUxtxf7Gn5TJzTmZ2ZGZHe3t7A0qWJKl11BXUEdFGLaSvycwfAGTmyszckJkvAN8ApjavTEmSWlM9V30H8C1gcWZ+qcv08V2avQN4sPHlSZLU2uo5R30kcArwQEQsqqZ9EjgpIiYDCTwCvK8pFUqS1MLquer7TiB6mPWTxpcjSZK68s5kkiQVzKCWJKlgBrUkSQXbqhueSGod9d74wxuESM3liFqSpIIZ1JIkFcygliSpYAa1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBTOoJUkqmEEtSVLBDGpJkgpmUEuSVDCDWpKkghnUkiQVzKCWJKlgBrUkSQUzqCVJKphBLUlSwQxqSZIKZlBLklQwg1qSpIIZ1JIkFcygliSpYAa1JEkFM6glSSqYQS1JUsH6DOqIeFlE3BERiyPioYg4v5o+NiJuj4gl1c+XNr9cSZJaSz0j6vXARzLzAGAa8IGIOBD4BDA3M/cD5lavJUlSA/UZ1Jm5IjPvrZ4/DSwG9gJmAVdWza4ETmhWkZIktaqtOkcdEZOAKcDdwO6ZuQJqYQ6M62WZ2RGxICIWrF69emDVSpLUYuoO6ogYDfwr8MHMfKre5TJzTmZ2ZGZHe3t7f2qUJKll1RXUEdFGLaSvycwfVJNXRsT4av54YFVzSpQkqXXVc9V3AN8CFmfml7rMugU4rXp+GnBz48uTJKm1jayjzZHAKcADEbGomvZJ4AvA9RFxJvD/gHc1p0RJklpXn0GdmXcC0cvsmY0tR5IkdeWdySRJKphBLUlSwQxqSZIKVs/FZJKGSGdnY9s1w1BuW2oFjqglSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBTOoJUkqmEEtSVLBDGpJkgpmUEuSVDCDWpKkghnUkiQVzKCWJKlgBrUkSQUbOdQFSFJJOjvLXp9ajyNqSZIKZlBLklQwg1qSpIIZ1JIkFcygliSpYAa1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWB9BnVEfDsiVkXEg12mdUbE8ohYVD3e1twyJUlqTfWMqK8Ajulh+pczc3L1+Eljy5IkSVBHUGfmfOCPg1CLJEnqZuQAlj0vIk4FFgAfycw/9dQoImYDswEmTpw4gM1J6k1n51BXIKlZ+nsx2deBfYHJwArgi701zMw5mdmRmR3t7e393JwkSa2pX0GdmSszc0NmvgB8A5ja2LIkSRL0M6gjYnyXl+8AHuytrSRJ6r8+z1FHxPeBGcBuEbEMuBCYERGTgQQeAd7XxBolSWpZfQZ1Zp7Uw+RvNaEWSZLUjXcmkySpYAa1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBTOoJUkqmEEtSVLBDGpJkgpmUEuSVDCDWpKkghnUkiQVzKCWJKlgBrUkSQUzqCVJKphBLUlSwUYOdQHSUOnsbGy7Zmxb6slQ/u5q8DmiliSpYAa1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBTOoJUkqmEEtSVLBDGpJkgrWZ1BHxLcjYlVEPNhl2tiIuD0illQ/X9rcMiVJak31jKivAI7pNu0TwNzM3A+YW72WJEkN1mdQZ+Z84I/dJs8CrqyeXwmc0OC6JEkSMLKfy+2emSsAMnNFRIzrrWFEzAZmA0ycOLGfm5OkbVNnZ2PbqfU0/WKyzJyTmR2Z2dHe3t7szUmSNKz0N6hXRsR4gOrnqsaVJEmSNupvUN8CnFY9Pw24uTHlSJKkrur5etb3gZ8Dr4yIZRFxJvAF4M0RsQR4c/VakiQ1WJ8Xk2XmSb3MmtngWiRJUjfemUySpIIZ1JIkFcygliSpYP294YkkaZjwpixlc0QtSVLBDGpJkgpmUEuSVDCDWpKkghnUkiQVzKCWJKlgBrUkSQUzqCVJKphBLUlSwQxqSZIKZlBLklQwg1qSpIIZ1JIkFcygliSpYAa1JEkFM6glSSrYyKEuQCpdZ+dQVyCplTmiliSpYAa1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMG55I0jDlzXqGB0fUkiQVzKCWJKlgBrUkSQUzqCVJKtiALiaLiEeAp4ENwPrM7GhEUZIkqaYRV32/ITMfb8B6JElSNx76liSpYAMN6gRui4iFETG7EQVJkqS/Guih7yMz87GIGAfcHhEPZ+b8rg2qAJ8NMHHixAFuTpKGp+F2c5J6+zPc+t0MAxpRZ+Zj1c9VwE3A1B7azMnMjszsaG9vH8jmJElqOf0O6ojYMSJ22vgceAvwYKMKkyRJAzv0vTtwU0RsXM/3MvN/N6QqSZIEDCCoM/O3wGENrEWSJHXj17MkSSqYQS1JUsEMakmSCtaIW4hK/daM71D6vUxJw4kjakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBTOoJUkqmEEtSVLBvOGJhh1veCJpOHFELUlSwQxqSZIKZlBLklQwg1qSpIIZ1JIkFcygliSpYAa1JEkFM6glSSrYNn3Dk625scW2cBOMemtsdLutbSupNfn/xNBwRC1JUsEMakmSCmZQS5JUMINakqSCGdSSJBXMoJYkqWAGtSRJBdumv0fdDM34jnKjteq2JbWuVv5/zxG1JEkFM6glSSqYQS1JUsEMakmSCmZQS5JUsAEFdUQcExG/ioilEfGJRhUlSZJq+h3UETECuBz4W+BA4KSIOLBRhUmSpIGNqKcCSzPzt5n5PHAtMKsxZUmSJIDIzP4tGPFO4JjMPKt6fQrwmsw8r1u72cDs6uUrgV/1sLrdgMf7Vci2z763plbuO7R2/+17a+re95dnZns9Cw7kzmTRw7TNUj8z5wBztriiiAWZ2TGAWrZZ9t2+t6JW7r99t+9bayCHvpcBL+vyegLw2ADWJ0mSuhlIUP8nsF9E7B0RLwHeA9zSmLIkSRIM4NB3Zq6PiPOAfwdGAN/OzIf6ubotHhof5ux7a2rlvkNr99++t6Z+973fF5NJkqTm885kkiQVzKCWJKlgQxLUEXFJRDwcEfdHxE0RsUsv7YbdLUoj4l0R8VBEvBARvV6qHxGPRMQDEbEoIhYMZo3NshV9H477fWxE3B4RS6qfL+2l3YZqny+KiG364sy+9mNEbB8R11Xz746ISYNfZfPU0f/TI2J1l/191lDU2WgR8e2IWBURD/YyPyLisup9uT8iDh/sGpuljr7PiIgnu+zzf6xrxZk56A/gLcDI6vnFwMU9tBkB/AbYB3gJcB9w4FDU2+C+H0Dtxi/zgI4ttHsE2G2o6x3svg/j/f4/gE9Uzz/R0+98NW/tUNfaoP72uR+Bc4H/VT1/D3DdUNc9yP0/HfjqUNfahL4fDRwOPNjL/LcBt1K7F8c04O6hrnkQ+z4D+NHWrndIRtSZeVtmrq9e/oLad7C7G5a3KM3MxZnZ093Zhr06+z4s9zu1PlxZPb8SOGEIaxkM9ezHru/JjcDMiOjpRkrbouH6e9ynzJwP/HELTWYBV2XNL4BdImL84FTXXHX0vV9KOEd9BrVPV93tBfy+y+tl1bRWkcBtEbGwug1rqxiu+333zFwBUP0c10u7URGxICJ+ERHbcpjXsx83tak+uD8J7Doo1TVfvb/Hf1cd/r0xIl7Ww/zhaLj+G6/X9Ii4LyJujYiD6llgILcQ3aKI+D/AHj3M+lRm3ly1+RSwHrimp1X0MG2b+C5ZPX2vw5GZ+VhEjANuj4iHq09rRWtA34flft+K1Uys9vs+wE8j4oHM/E1jKhxU9ezHbXZf16Gevv0b8P3MfC4i3k/t6MIbm17Z0BvO+70v91K7x/faiHgb8ENgv74WalpQZ+abtjQ/Ik4DjgNmZnXwvptt9halffW9znU8Vv1cFRE3UTuUVnxQN6Dvw3K/R8TKiBifmSuqw3yrelnHxv3+24iYB0yhdq5zW1PPftzYZllEjATG0ITDhkOkz/5n5pouL79B7XqdVrDN/hsfqMx8qsvzn0TE1yJit8zc4h8qGaqrvo8BPg4cn5l/7qVZy96iNCJ2jIidNj6ndvFdj1cRDkPDdb/fApxWPT8N2OzoQkS8NCK2r57vBhwJ/HLQKmysevZj1/fkncBPe/nQvi3qs//dzsseDywexPqG0i3AqdXV39OAJzeeFhruImKPjddhRMRUahm8ZstLMWRXfS+ldo5iUfXYeOXnnsBPurR7G/BraiOKTw1FrU3o+zuofaJ8DlgJ/Hv3vlO7UvS+6vFQK/V9GO/3XYG5wJLq59hqegfwzer5a4EHqv3+AHDmUNc9wD5vth+Bi6h9QAcYBdxQ/X9wD7DPUNc8yP3/fPXv+z7gDuBVQ11zg/r9fWAFsK76934m8H7g/dX8AC6v3pcH2MK3X7a1Rx19P6/LPv8F8Np61ustRCVJKlgJV31LkqReGNSSJBXMoJYkqWAGtSRJBTOoJUkqmEEtSVLBDGpJkgr2/wH01ScgGtUbRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x13dcb0c5240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#在训练集上观察预测残差的分布，看是否符合模型假设：噪声为0均值的高斯噪声\n",
    "f, ax = plt.subplots(figsize=(7, 5)) \n",
    "f.tight_layout() \n",
    "ax.hist(y_train - y_train_pred_lr,bins=40, label='Residuals Linear', color='b', alpha=.5);#（真值-预测值）并构建直方图\n",
    "ax.set_title(\"Histogram of Residuals\") \n",
    "ax.legend(loc='best');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "通过直方图可以看出，残差大致呈正态分布，但是略微向右偏移，并且基本可以确定左侧多余的点为噪声点（因为其残差绝对值较大）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAEYCAYAAABiECzgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXl8VPXV/99nkgmZsIUliglr0YKCsqWKYktRK64YxaK4tNantfbRVm3FYutP0McFxaXY2seitRvUoiJRShFFUJ+KqAmggEJdECGIJELYMoRJ5vv7YzJxMpk7c2fLzWTO+/XiRWbmzr3nziSfe+75nkWMMSiKoihtj8tpAxRFUbIVFWBFURSHUAFWFEVxCBVgRVEUh1ABVhRFcQgVYEVRFIdQAVYURXEIFWBFURSHUAFWFEVxiFynDYiH3r17m4EDBzpthqIoSlQqKytrjDFFsbbLKAEeOHAgFRUVTpuhKIoSFRHZamc7DUEoiqI4hAqwoiiKQ6gAK4qiOIQKsKIoikOoACuKojiECrCiKIpDqAAriqI4hAqwoihKBBobG2loaEjrMVSAFUVRwqipqeHss89m+vTpaT2OCrCiKEoYjz/+OK+//jpDhw5N63Ekk6Yil5aWGi1FVhQlHRhj2LlzJ0cddRQNDQ385z//4bjjjktoXyJSaYwpjbWdesCKomQ9dXV1fO9732P06NFUV1eTm5ubsPjGgwqwoihZzYcffsjYsWOZP38+P/nJT+jVq1ebHTujuqEpiqKkkvLycr7//e+Tm5vL0qVLmThxYpseXwVYUZSsxBjDk08+yZAhQ3jmmWcYMGBAm9ugAqwoSlbxxRdf4PP56Nu3L/PmzaNTp0506tTJEVs0BqwoStbwxhtvMGrUKK688kqMMXTr1s0x8QUVYEVRsgBjDHPmzOHb3/42BQUF/OY3v0FEnDZLBVhRlI7NgQMHmDp1KjfeeCPnnHMOFRUVjBgxwmmzABVgRVE6OMYYNmzYwKxZs1i0aBGFhYVOm9SMLsIpitIhWbJkCRMmTKBr165UVlY6Guu1Qj1gRVE6FD6fj5tuuonzzjuPhx9+GKBdii+oB6woSgdix44dTJkyhTfeeIOf/exnTJs2zWmToqICrChKh2D16tWUlZVx4MABnnrqKS699FKnTYqJhiAURekQFBUVMWjQIN56662MEF9QAVYUJYPZu3cvDz/8MMYYBg8ezKpVqxg2bJjTZtnGMQEWkX4islJEPhCRjSJyg1O2KIqSeWzYsIFvfOMbTJs2jbVr1wK0i+KKeHDSA24AfmGMORYYC1wnIulvwKkoSsYzf/58TjrpJPbv38/KlSsZPXq00yYlhGMCbIz53Bizpunn/cAHQIlT9iiKkhncdtttXHHFFZSWlrJ27Vq++c1vOm1SwrSLLAgRGQiMAt5y1hJFUdo748ePp76+nnvuuQe32+20OUnh+Ew4EekCvAbcbYx5LsLr1wDXAPTv33/M1q1b29hCRVGcZvny5WzcuJEbbsiMpaKMmAknIm5gITA/kvgCGGPmGmNKjTGlRUVFbWugoiiO4vf7ufvuuznzzDN58sknqa+vd9qklOJkFoQAfwQ+MMY85JQdiqK0T/bs2cMFF1zAbbfdxtSpU1m1alW7LSlOFCdjwOOAK4H1IrKu6blfGWP+5aBNiqK0Aw4fPszJJ5/MJ598wm9/+1uuu+66jEsxs4NjAmyM+TfQ8T5RRVGSJi8vj2nTpnHcccdx8sknO21O2tBKOEVR2gWHDh3immuu4fnnnwfgv/7rvzq0+IIKsKIo7YAtW7Ywbtw4Hn/8cTZu3Oi0OW1Gu8gDVhQle1m6dCmXX345fr+f559/nkmTJjltUpuhAqxkBeVrq5i9bDM7ar0UF3qYNnEIZaO08NJpKisrOffccznhhBNYuHAhgwcPdtqkNkUFWOnwlK+t4tbn1uP1NQJQVevl1ufWA6gIO4Tf78flcjF69GieeOIJLr30UgoKCpw2q83RGLDS4Zm9bHOz+Abx+hr5xdPvMmj6EsbNWkH52iqHrMs+3nnnHU444QTef/99RISrr746K8UXVICVLGBHrTfi843GYPjKI1YRTi/GGP7whz9w6qmncuDAgQ5X1ZYIKsBKh6e40BNzG6+vkdnLNreBNdlJXV0dV111Fddeey0TJkygsrKSUaNGOW2W46gAKx2eaROH4HHnxNzOylNWkmfOnDn87W9/Y8aMGSxZsoRevXo5bVK7QBfhlA5PcKEtmAXhEqExQhdAO55yOJpdEZ0DBw7QpUsXfv7zn3PKKacwfvx4p01qV6gHrGQFZaNKeGP6aWyZdS4PThnRyiP2uHOYNnFIXPsMZldU1Xo1lhxGQ0MDt956KyNGjGD37t106tRJxTcCKsBK1lE2qoR7LzqekkIPApQUerj3ouPj9lytsiuyPZa8a9cuJk6cyKxZszjjjDOyNsPBDhqCULKSslElSYcKrGLGqY4lZ1KY48033+S73/0uX375JX/605+46qqrnDapXaMCrCgJUlzooSqC2CYSS7Yi04pIZsyYQX5+Pm+++SYjR4502px2j4YgFCUC5WurGDdrRdRCjUjZFYnEkqORCWGOAwcOUFNTA8C8efOoqKhQ8bWJesCKEoZdrzM8uyId4YG2CnMkyqZNm5g8eTLFxcW89NJLHHHEEU6blFGoACtKGNG8znBxTUUsGazjvG0R5kiUZ599lh/84Afk5+czZ86cDjmxIt2oACtKGG3tdUbzuKdNHNLiNWgZ5kjVAp3d/ZSvreK+f23k/ecfY987i/j68aNZvqScfv36JXLqWY8KsKKE0dZeZzSP+43ppzVvU1XrJUek+bWKrbtZWFmV9AKd3ZBLcLsD+/dy8D+r6Dr6PGTiNVTWuFD9TQxdhFOUMNpicS1I+dqqiGIPASEcN2tFC5uCFXxVtV7mrf4sJQt0dhf6bn/ieeoOHSInvwvFV82h53eu5ZDf1a4WBDMN9YAVJYy2WFyDrzzKaAS90U65rlYiGe0942ataPaYG42hJMo5xAq5GGN46KGHWP+/t9B93FQKx03Fld8l5vuV2KgAK0oEUrW4Fo1InmckvL5G2+ILgVHjQa861GO2Ck9EC7ns27ePq6++moULF9Jz2DfpXHpBxO2UxNAQhKI4RLo8x9ZthgJ4fY3cuGBdq7xmq5DLpccIJ554IuXl5TzwwAM88de/07lL11bbpSM0ky2oB6woSZBMFoKV5xmJHgVuDvn8cXnCVoR7w1Yhl6Pde/iNz8crr7zS3EhHRDKmLDoTEBOhLV97pbS01FRUVDhthqIArbMHIOAR2m3sU762ipsWrLP0WIO4c4TZF48A4MYF65IxuQUlhZ7mLIsghw8f5tlnn+Wyyy4DAl3NcnNzm+1V8bWHiFQaY0pjbachCEVJkGTLhMtGlXD52P7EKl/onJeblph0eAhk+/btjB8/nssvv5y33noLoIX4auvN1KMCrCgJkoqCjbvKjufhS0aSE6WKbK/X1/xztO3iJXTxbMWKFYwePZoNGzZw832P8fOVB1v0wciEnhSZiAqwoiSI1ep/vFkBZaNK8EcJBRYXepqbA0Wa5JEIoYtnc+bM4Tvf+Q69e/fmnj8v5p/7B7TydK1i1ZqClhwqwIqSIKks2LASbQEmDC2KKoKJYbipKSPi84bOTJkyhbfffpt/fNh6oc/ra7T0vIN22+kep7RGBVhREsRqsgYQtxhZDQ41ELHiLVn2Vn3E/g2BYo3na/tyyS0P0KVLF0uPttGYVrHq4MVG48OJo2loipIE4YtjiTZQD00Fq6r1Iljn89qh0OOmNiR2HMqB9a+w+6VHcXXuQeehp+Ilr7nTm1VqXLg9AkweEzj3cbNW2O4ep7REPWBFSSFWi1W/ePrdmB5xcHBoSaEnKfG9Ymx/1s04k5KwsIZp8PHlst/x5b8eJq94KEdd+SCSmwd8Vb48sFfkUEi4PQZYuakaaP89i9sz6gErSgJY5cRGu4UHex5xMsJV6HFzV9nxzccKYhob2Pn36Rz+fDPdxl5M4TevRFwtQx5Vtd64jh3ctj33LG7vOOoBi8iTIrJLRDY4aYeixEP52iqmPfNui5jntGfepXxtlS3RiZW+lahwedw5zJw0rNnG0Jit5ORSMOQUii66jR7jr2olvkHi8byDdqZiMTJbF/GcDkH8GTjLYRsUJS5mvrARn7+lVPn8hpkvbLRcTAsnmqdpdx/hhFbgzV62Gb/xU/vvv+P9NFA91/2kyXQ+ZixXjO3fKjwRL6ECa7UYaTf+m82LeFFDECLiAi42xjydjoMbY14XkYHp2LeipAurxa1ar69VXwVXUzvIcKJ5ueELcqEtJQ/WN0Q8flBQx81awY5aLw3efdQsfpBDWyrpWnoBnoGBIZkGmkMUwZaVdskRwW9MxDLk8PMOevh2RDieEVAdjagCbIzxi8j1QFoEWFEyDTteWWhmhFW/iFi35+H7CApbYYEbt0taeOAedw4De3ma+0rU7/yI6kX30HhwNz3P/G+6jDy7eduSpqKOoLjHQ6Mx5IhQVetl5gsbuWPxRmrrfBQXepgwtCjh6RzZXORhJwTxsojcLCL9RKRn8F/aLWtCRK4RkQoRqaiurm6rwypKRKLFbnsUuFs9l+rb8z11PpDAYltwf5PHlLDq490Y4PCuLeycdzMYQ5/L7qPrqHOah2W6c4Rd+7zcuGBdwkUdQW++1utjT52vOWQwP8HpHOGx6lCyYRHPThbE1U3/XxfynAG+lnpzWmOMmQvMhUA3tLY4pqJYEc0rm3H+sIjPJ9NIJ9Ltua/R0LlTLutmnAkEQgnBPwx30UAKx11GlxETySno3vwel0Cj3+BP01+Q1W5jebGzl22O+F6BrOgzHFOAjTGD2sIQRckErFKuCj3utMQro+XYBkMJn366hS+XPkKvs3+Gu7AP3U+e0mr7dAlvLGJ5sVbnZ4hvsGimEjMEISIFInKbiMxtenyMiJyXioOLyFPAm8AQEdkuIv+Viv0qSrqIlKHgdgkipCWFykrACgvc3Prcej6seJXP/3wDvl2f0LivfYXo3C6J6cVanV+yWRqZQsyG7CKyAKgEvmeMGS4iHuBNY8zItjAwFG3IrrQFsRqPh77e3ePm4OEGfI0t/44KPW5mThoWlxd3W/l6nnprW/Ni19ST+lE6oGfERbw8l+HTl//MvjefJu/IwfQuuxV3YR/LfUcrTU4XPQrcrL39zKjbJNvUvr2Syobsg40x9wM+AGOMF2L2kFaUjMROTmqwZHjLrHPp3Cm3lfhCYJEqnlzW28rXM2/1Z82LXI3GMG/1Z1Rs3R1xEW/b60+z782n6XLCmfS5YnZU8RVg5qRhuHMS/7MNNkMLdkUr9LjpUeCOKgS1dbEFP9lFykzHziLc4Sav1wCIyGCgPq1WKYpDxJuTGm2RKZ5c1qfe2mb5/F1lXwnSwne2MnvZZrqMOo+cbkfS+dhvxty3AX7x9Lt8raiAD3cdjLm91U4+nXVuxJes8ontZjG0xQTq9oodD3gG8CLQT0TmA68At6TVKkVxiHgbyyS6yBSOVaP14PPGGK6ZfjeXn38a23btxpWXb0t8Q/eTsPjSsil8eKw7lX2Rs42YAmyMeRm4CLgKeAooNca8ml6zFMUZ4p1yEats2K4XaNXwPEeEgwcPcuWVV/L4fbfh6tIL/KntDRwLjzunRVP4YGjmpgXrGDh9CbOXbWbymJKsDSMkg91uaOOBUwnczbiBRWmzSFEcZNrEIXFVrgVF5o7FGwNFEiGEvy/a4t7Uk/oxb/VnrfZ/Zt9Gxo4dy8aNGyn85hV0O3kKgQ4BbUNJk52RQjNBn72q1svCyioV3QSIKcAi8nvgaALeL8CPReQMY8x1Ud6mKBk5xjy8p4EduyOJcHgWxG3l65m/+rMWohVaqhvszxCeBbH6dzfx+eef8+KLL3LHmtwUjyWyJjwT4aYF66Juny29G1KNHQ94PDDcNOWrichfgPVptUrJeBKdDNEeiHdRKFIqVX2Dv8XroeIbxOtr5I7FG5uPdVfZ8dxVdjwNDQ0cOnSILl268HjePdS+/jE/XuGjsCCp0wJg3OCevPnJ7oiFGdGa7VgVoISSDb0bUo2de5nNQP+Qx/2A99JjjtJRyKYx5jNf2Bj1XK3KbSHQ2yE0VW3nzp2cfvrpXHbZZSxas52H3txNDd2a+0C4kkgAdQkMKuoSMd7szhEenDKCLbPO5Y3ppwEt59pNGFoUs0Vmqns3ZEOPYEsPWEQWEwjzdAc+EJG3mx6fBKxqG/OUTCVbxtSUr62yLHAInqudfggA/+8Pz/LB/Dsx9Qf56e3388BL/2kl7MmUFPvNVyGOcDrn5UYNlyysrGLymBJWbqqOOLMu1VkPmXwHFQ/RQhAPtJkVSocjW8bUxJpsUb62yrIncJDte+r48S/vZNfyJ8gtPJIjptzJsvrBeA+k/mJlZcfepotItHDJyk3Vzd5xMvF9O+/Nlh7BlgJsjHkt9LGIdIu2vaKEEm82QaYSzbsNpm5FE18AqT9Azaqn8Rx9Ir3PvQlXp854fY3NjdhTidU+gxfGaOGS0HNNtHjCrmebLXdQdprxXCMiXxCI+1YQ6AuhDRmUqGRLiamVRy8C8yL0yA3FV7sT42+E/K4cdeVDFF34a1ydOje/3mhMq7hrMj0A3K5AZkW0ooloAucSSToea3dtIN587EzFjkc7DRhmjKlJtzFKxyLVJaZtkdYWqSFOMEUsEpE8fYBYjuvBD17ny6WP0P3kKXQ/eQq53Y9otU2OCJPHlLQICeTluvD7TauZdLYQKB3Qk9IBPVuMOwoVwGjZDvFMdrbCrmebLXdQdrIgPgbq0m2IokSjLQY3WjXEua3cOusy3NO3qmgLYhp97F4+l5oX7ifviEF0Hn6a5baNxrDg7W0tQgL1DX78BPKM48XXaJpjqMEKvnBRtZPtAIlntNj1bLPlDsqOB3wrsEpE3iKkCY8x5mdps0pRwmiLRZlYDXGsCPX0B01fYrldw/4aap6/j/qqD+haegE9vv0DJMf6TzBHJKKn2+gPTMTY6/XFNUYevvI0rT7PlZuqufei41vcaaRyZls8nm02NOmxI8B/AFYQKL7wx9hWUdJCWyzKxGqIY4doguU/8CW+L7fTe9ItdD72W1H343HnRI0fh05Ljoegpxnt8wwXvmS7nYWSSKVhR8aOADcYY36edksUJQrR0tpSFRu2ErRYYYVQwj08YwwN29ZTMGgEHDWEkmv/iKtT5JK28Eq0WJOL4xXfUE8znjTBVMdjs8GztYudGPDKpkyIo5yYiqxkJ+FVUJFik1ZduhKNDU89qV9cz0eibFQJk8eUBMS0/iDVi+5mx1O/Yv/WjQCW4gtfCerB+gbuWLwx5X0fQmOo8bSQzJZ4rBPY8YAva/r/1pDn2mwqspJ9RMoVDa3ECvV0UxkbtmqIYxX/jeR5AyysrML7xSdUl99Dw95d9Dj9R3QqOdaWDQbSMjoo0tDQTrmu5s+uR4GbGedbj1BSrzU9xJwJ157QmXDZgVXMsaTQ01yJFWRglEUvqwkO8WAV3rCaZeYS2LXuFb5c+ltc+Z3pfcF08vsel7QdwSGViXrFwSBKcaGHCUOLAheJDjaHrT1hdyacnXaU34v0vDHmr4kYpiixiGfBLRVxWyuiVW1Zed4Axt9IXvHXKZp0CzmdeyRtBwTO/eFLRkbMObZDaF+HSH2HO2KZbyZgJwTxjZCf84HTgTWACrCSFuJZIEpF5oIV0cIb4ReDhn27OLxrCwVHn0SX4afTediElDZOLy70tMogiNVjIl46WplvJhBTgI0xPw19LCLdgb+lzSIl64ln1b3EQqxLUlCyaiVI4Slg3i1rqFn8ALhc5P94BC53fkrFV6D53ENjsZHCIMnQ0cp8M4FEfkvqgGNSbYiiBIln1T2dAyGjCVKjMRjjp/aNp9j19AxyOvegz2X34XLnJ33ccKx83ODn1KMg/qq4cDpimW8mYCcGHOwLDAHBPg54Op1GKYrdVfd0JvZb9XmAQJy3+rm78H78Dp2HTWD4lF+w1+fC60tPrZJV74WyUSXMXra51Tw6O0SbgKG0DXZiwKF9gRuArcaY7WmyJ+vJxDlqTpOuFKlQcQ8Pc4grB3fRADxfK6Vm2aM8v24Hdyze2EqA3S7hkhP7tUifSySTIdoiWbTYrQD57tYXBrdLmP3dEfq75TB2YsCvxdpGSQ3ZMgXASUIvcIUFbowJNCO3utgFxT2YGrf/3ZfI692fTiVD6TH+KkoKPTy/bkdET1kELjmxH6UDerJyU3Xz8z0K3Al5rFZCayXqwYnGgZFJYZ558kkiSgqwE4K4CLgPOILA1yaAMcZ0S7NtWUe2TAFoa4KiGz5KJ1QEI13sytdWNU869vvq2f3yYxxc/zKdh59Op5KhzXHTSN8bBFpSLnhnGwve3tbcVKeq1pvwXLfQmHTohaS7x407R/A1fhUtDq0SjGRbaGc0xTnshCDuB843xnyQbmOynWyZAtCWhN9VREvaCr3Yla+tYtqz7+JrNPhqd1JTfi+Hv/iY7idfQvdTL2tRORZtZHuoKAZJsJUvVbVeBt/6LxqNaXEhqfX6cLuEHgVuaut8UasEQ9HfK+exI8BfqPi2DdkyR60tiSVC4VTVehk3awUH6xsC4vvlNnb+7WYAii6eQcHgQFp8QcgQy0TjuvEQFNtg6lu4hvv8hoK8XNbefmbzc9EuDKC/V+0BOwJcISILgHJa9gN+Lm1WZSnZMgWgLUnEywsV09wexXQefjpdSyfhLuzTYptB05c0l/ZGGmQZLz0K3Bzy+Vt8/+HTh6MRfq7RLgz6e9U+sCPA3Qjk/p4Z8pwBVIBTTLIpVZmaQZEquyPtJxHvtLFuL7tfmUuPb19Nbtde9DzjmojbBbuvLaysikt83TkChhbN1j3uHGacPwzAVjP0SIR7tFZpdLEa7yhthzbj6SBYNYdp7w1WUmW31X4mjylp1Xgm6FW6pHU8tn7HZqrLZ9FYV0vRBdMpOOYkW8e32xw9mJkA9i60Vo2JwrH6zDL1opzp2G3GowLcQYing1h7IlV2R9tPcEEqKEKRuoEZYziwbim7l88lp2svispupVOfo+M6l1hTLBL5LuyUG9sZHposKuTxkbJuaOlERM4C5gA5wBPGmFlO2pPJZGoGRarsjnfETrig7V/zT/Ys/wP5XxtD7/NuJsfTNa7jhwp9eLobJB5zDS8GCXraoftvNIaFlVWUDuiZFlHU/PT04ZgAi0gO8CjwHWA78I6IvGCMed8pmzKZTM2gSJXdVvvxuF2Mm7Wihfcbup0xBhGhy7AJYPx0HXN+i0Y6dhfBgh5haA6xXY8x1rZ2ZrSlM19c89PTh6UAi0jUOXDGmIeSPPaJwEfGmE+ajvcP4AJABTgB0pFB0Ra3namye9rEIUx75t1WU4TrfH7qmsQqvBdu3X/eZP+af3LExTNw5XehW+kFrfZrgCvG9o/YQzeU4Ij24OdjtzzarncZ+l1YXRB21HrT8p1l6t1VJhCtG1rXpn+lwE+AkqZ/1xJoyJMsJUDoHPDtTc+1oGkeXYWIVFRXV4e/rDSR6rldQWFIxay1trC7bFQJXfLt3dAZfyN7Xv0T1Yvuxl9fh//Qwajbz1v9WcwG74l+PtG8yyDh34UV3T3utHxnVncj7f3uKhOw/I01xtwBICIvAaONMfubHs8EnknBsSP9Rrf6/TLGzAXmQmARLgXH7bAk2pQmktfUlredqWqmU2ujv0LjwT1UvzCb+s/eo8vIs+l5+jVIbux2jnYyHBL5fOx4l3aKSTzuHERIy3em+enpw47L0B84HPL4MDAwBcfeDoSOm+0L7EjBfpU4sLoFtvqDT/dtZ7Rb6Fi313byZmv++RCHd2yi17k30WX46UAgi6Brfq6tYZjBFo7RwgB2zyeazcWFnhY9LKJR6HEzc9Iwy8q3ZL+zdLb8zHbsCPDfgLdFZBEBD/VCUjOO6B3gGBEZBFQBl/LVBGaljbDydK3yWtN52xktHgq0eu2mBeuo2Lq7Of3KqvDAGAP+BiTHTc8zfoxpPEzeEV8N9W40hoOHG3C7pFUMOZxgBoKdz8fqfCq27m5uT5lII51w6hv8zcdO10KsTkVODzEnYhhj7gZ+AOwBaoEfGGPuSfbAxpgG4HpgGfAB8LQxZmOy+1Xiw8o7ajQmbZMmrAj00418Cx3pQmGA+as/a45xBuPJofgPe6l54X6+/NccjDG4e/VtIb5BfI2GLvm5toZ5GqxDEgfrG5rtsbq4zV/9WXOcttbrw9do6JyX0yIGvnJTte0eFsHPKJ3TQZT0YHckUQGwzxgzB9je5LUmjTHmX8aYrxtjBjcJvdLGWHlHQSFI1aJeLMrXVln2yN1R67W8UBhosWBVNqqkeR6c78tt7Pzrz6nb/AbuogExbait8/HglBGtRCwa4XJd6/U1L3xZhQ4iSXfd4UYevmQkb0w/jbJRJXGXTwfzndvyO1OSx04/4BkEMiGGAH8C3MA8YFx6TVPagmgLLMnedpavrWLmCxubY6vRehCEimg4wYuElSiFi/O0iUO47q7f8/k/H0Zy8zhiyp10HjQyZhvIwgI3ZaNKqNi6m6fe2mZr4S3SFkGP1G55cnA/oa0w42nCA199RhoqyCzsxIAvBEYRGEWPMWaHiMRXJqS0W9K1wFK+tqpVXu6eOh+/eObdFjHQ4PGiLRQFb6FvtFhk6u5pmcUwfmAB+5Y/Spc+g+h+7i3079+Pg/UNMRfZDhxq4Lby9SysrGohnB53Dp1yXbYW6YJEy9eN9h4IfBfxvDd0arKSWdgR4MPGGCMiBkBEOqfZJqWNSYfXNHvZ5ogLWo1+06KoIbgw1d3jjihwhR53s22/eu496iIMvQyGbWtqaujVqxc9evRg1b//j03eLvxmxRbbYujzm4ier9fXSL47vgHisbz2aO+JN2vBoCXBmYqd36qnReQPQKGI/AhYDjyRXrOUTCceEfH6GhEh4gLSzEnDQraLPHG4ts7Ha6+9xvDhw3nooUCB5scNPbl98eaYxQvhWIUMaut8tseoBT3SSItiAHk5rfcUulgWb9ZCiRZEZCx2siAeAJ4FFhIazK1HAAAcG0lEQVSIA99ujHkk3YYpmU28IlJb52u1gDR5TGDk+qDpSxg3awWFBa0LJowxmPde4PTTT6ewsJCzzjoLiH8SRhCrLIjiQo9tIQ96pMFFsR5hdh9uNLhzhEKPO+JimZVwR0KzHDIbO4tw9xljfgm8HOE5JcNJV7+HaROHWMZsI1Fc6GnVzCY8h9btkhY5s/76g+xZOocDm1dx8cUX8+STT9K1a2B5IpHiA6v+waHDN+2EFEI90rJRgYtIeIZH8By2zDq31futOqCVNDUTCo+fa/ghc7ETA/4OEC62Z0d4Tskw0tlmsGxUCb9etJ6Dh2N7oZG8uEgerM9vKPS46dwpN1DEULeDnVsqeeihh7jxxht5ft0OZi97hx21XlxxZCAECXqhpQN6RrwoVWzdHbMpT6RzsboY1Hp9lK+tivhZazZDdhCtG9pPgP8GBovIeyEvdQVWpdswJf2ku9/D3Rce3zxZOBSP20W+O6fFBN/w41mJ1l6vj3mTixk+fHhguxlTKC4ubnUxiVd8S5o8cLAWv5WbIjeDCpYnW51LtBJpbemY3UTzgP8OLAXuBaaHPL/fGLM7rVYpbUK62wwGheWOxRtb3IIHFtOEhy8ZGVF8ytdWRfRgTYOP+jf+xAn3L+b111/n1FNPpbi4GEg85gv246hWn4vfmIihhCDRwjHa0jG7sVyEM8bsNcZ8SmBixW5jzFZjzFbAJyL2BmUp7Zq2aDNYNqqEgrzW1/nwlotBgp5suPg27Kvmi6ems3P1C3T9xoVMe62uueQ3WtVZOB53DleM7Z9QtViin1fZqJJWC3F236t0bOykof0vcCDk8cGm55QMp616B8TjaUfsB/HpOj7/8w0crvmM3mW30mPC1Xy+P1Dye1v5+hYNe8IRoVVmRaKLWMl8XjPOH6Z9GpRW2FmEExMyudMY4xcRR2fJKfawM+oGolfBpSJLwm6XLqt+EL4vt5FTUEjRhb/C3atv8/NeX2PMkmFjaB6EmeyiYzJVg9rSUYlEzKnIIvIc8Cpfeb3/DUwwxpSl17TW6FRk+6Ri3Huq9hEeA4bArVf3AneLhbjQNK/GQwfw1Wwlv++wQK5vw2Fc7k62jhlK6CTiTJ0crWQedqci2wlBXAucQqBn73bgJOCa5MxT0o2dUTfp3kdQwCN5tX4CvSFCR+cExfHwF5+w8y83Uv3c3fgPexERenbrEvEY0dpHht/i62wzpb1hpxJulzHmUmPMEcaYI40xlxljdrWFcUripEJskt1HPJkJwSbwB957mZ3zbsY0NnDE5Ntx5XmaJz5EiqFOPalfxKqxQo+7laeus82U9ka0POBbjDH3i8hviTyr7WdptUxJilRMR0h2H/GIvfE3suul33Pg3WXkDxhB70m3kFPQvbkfRLQYqlXhRDhtPdusLaZKK5lNtMW0D5r+16BrBpIKsUl2H3ZmtAURVw4FeTmcMPmHfD54En5xkSPC5DElMQsk7FaNteVCWDqrDJWOQ8xFuPaELsLFRyo8sNvK1zdnGuSIMPWkfs0z2OwcP9ZcM+/HFeR060334sFcNLqY59bsiLjoB5mVQaALftmN3UW4aCGIxURpym+MmZSgbUobkYqJFqHNyRuNYWFlFaUDesbtcYaLkfE3sveNf7B31T/oOfxU7v3rU5aLfncs3sghn9/Sm2yPt/q64KfYIVoI4oGm/y8C+hAYQwQwFfg0jTYp7YRogmhX+IIXgVCPsNG7j5rFD3BoyxqKRp/Jp/+3iIKCAsty3UhZFKHZGO3xVj+dE4qVjoOlABtjXgMQkf8xxnwr5KXFIvJ62i1THMfKW9tT52Pg9CUt5pZFEr5QgQ6OX6/78nO+eOpWGg/uoc85P+X3/3MLBQUFCc1B21HrTXtDoURp6wU/JTOxU9FWJCJfM8Z8AtA0EbkovWZlLu3xdjhRYi2ihYtlqPCFz4Sr9fpwAb2P6MPekmP52oRLuPOHF7QIU0QSXwHLcUUuEduDOtsarXxT7GBHgG8CXhWRT5oeDwR+nDaLMpiOtvIdb1N1+Er4Zr6wsVl8/b5D7P2/+XQ7+buIpxcH3n/N8n3hGGDmpGERF/OilSC3h1t97emrxCKmABtjXhSRY4ChTU9tMsbUp9eszKS93g7HQ7gHX+B2RRyEaUVwQnHQY/Xt2UH1onvwVW8lr89gaj3fZtysFa28QStvO7RPb9CuWM3W9VZfyRTsjCQqAH4ODDDG/EhEjhGRIcaYf6bfvMwikZXvtg5ZRDue1Rggl0CEAccRCa0MrvvwLWqWPISIiyO+OxPP18Y07zf8ziBWzDTUmxw0fYnl8Uv0Vl/JIOyEIP4EVAInNz3eDjwDqACHEe/Kd1uFLIKiW1XrjbpwZjUGKB5qmzIW/JtXUl3+IHl9jqao7FZyux/ZYrvwO4N4YqbRvGXNsVUyCTvNeAYbY+4HfADGGC/YntCdVcTbLzYVDXNiERT5oGBZLZxBahaughebO6+/ksKxk+lz+f2txDdI+PHKRpXwxvTT2DLrXN6YfprlRait+hgrSrqx4wEfFhEPTX+7IjIY0BhwBOJd+W6LZP1IDc6tjhdP6XAkzM7N5Lz7Og03f4urTh9BYc85URfxigs9CYVgNMNA6SjYEeAZwItAPxGZD4wDrkqnUZlMPCvf6U7Wt2pwbnW8CUOLYk79jYQxhsb1S9n58lxy+/dnx44d9O/fn7JRJcx8YWPEFDJpOl6iIRjNMFA6AlFDECIiwCYC1XBXAU8BpcaYV9NuWRaQ7lvpmS9sjLlN6PGspv5Gw3/YS83iB6ha+nvyBo6m9xUPsebLwDmVr63i4OGGiO+7fGx/Vm6qTnsIRlHaM1E9YGOMEZFyY8wYwHrpWUmIdN5Kl6+tiuh5huJxu3AJ3LhgXcx8X487h3y3q5VHXbN4Nt6PKyj81vfoNvZivqh3NXuxs5dtbjWSHqBHgZu7yo63zGZwuohCUdoKOyGI1SLyDWPMO2m3JgtJ1a10aCy1sMBtK/TgtZnfmyPS3JEsWN1mjEFEKDz1crqOmYRn4MiQ/TY22xKJYKaE9ktQsh07WRATCIjwxyLynoisF5H30m1YplO+topxs1YwaPoSxs1a0TxCPV3HCmY6GCI3r0kUjzuHB6eMaL5IGH8Du1c8we6XHwMg78jBLcQ3SNCjj0Tw+UghmGBsWFGyATsCfDbwNeA04HzgvKb/E0ZEvisiG0XELyIxe2ZmGuGCGFxcSpcIxzP6Jx5KCj0txvrc/cwqts//FfvfKUdEMMbag3aJMGFoUdQYd9moEiaPKWmR02iAhZVVab1gKUp7wVKARSRfRG4EpgFnAVXGmK3Bf0kedwOBhb0O2VWtLfJ7Q0lHzDRY1BAU39dff521j/yYw198RO/zb6bnd65FxPr6HewdPHlMCSWFHoTWgg6Bhb9oucmK0pGJFgP+C4Hii/8j4AUfB9yQioMaYz4AkCgTbTOZRPN7Ey1LTjZ/NxKhtu7du5dJkybh9nSlxyX/Q17RwBbbFnrc7D/U0Ko/g9fXyMpN1VGr07RxuZLNRBPg44wxxwOIyB+Bt9vGpJaIyDXANQD9+/d3woS4hdHu4lJ4v9yDhxuaswaCYYuKrbtZuak66rETzd+NRnGhB6/XS35+Pt27d+eWB5/kqY+E/SavxXZulyBi3ZkslpDqQpySzUSLATev5BhjIidzRkFElovIhgj/LohnP8aYucaYUmNMaVFR2y/OJBLPtZPfG77fWq+vVcqW19fI/NWftTj2tGfeZdSdL7VY3EskfzcaAlxyNIwcOZK5c+dSvraKv31a0Ep883IEn99EXfSLJaRaVqxkM9E84BEisq/pZwE8TY+FQIpwt2g7NsackSIbHSWRFpN28nvtLpyF+5WhgldV62Xas+9GzLVNFAFG+jZw6/fvpGvXrgwdOpRfWdh6OMZx7QiplhUr2Uy0kUQ5Vq9lE4nGKGPl96YqxpkK8e3RlDfs8jdQs+KPlFcu5tiRJ7J8ySKKi4vZsTT+Gpx42kJqWbGSrdhJQ0s5InKhiGwn0OJyiYgsc8IOO8TKZU31fp1gxvnD8LhzqNv2Pvsr/0nXb5ThP/v/8fYXAXGP19bwDApFUSLjiAAbYxYZY/oaYzoZY440xkx0wg472I3nxlt0EWm/kWiLPJG7n3kDr6+R/AEncNTVv6PnaT/kkF+aU8EmDC1qZYeVXQIav1UUm9gpRc5qYsUoYzVVt8qgCN9veBYEBIR+8piS5iyISNskgzF+fJWLqHh1Hn2m3kOnkmPJKxrQ/HpVrbfV9GMIiOwpg3uy5rO9LWLDQqDJjnq+imIPMVFma7U3SktLTUVFhdNmtGDcrBWW0xkmDC1i/urPWoiXx53TqhghiJ10t/K1VdyxeGPc5cZul+AHGoODMg8dYPe/Hubgh29ROHw8Xb9zPa48+6GGYIxXF88UpTUiUmmMiVnlqx5wklgtplXVeluJL0TPoAj1jINifNOCda3E7VAcQzIhUCgxc9IwIOBxf/qfjex+fha+fbt45JFHeOLLr7P3UHyZhjtqvbp4pihJ4kgMuCNhtUCVI9JKfIPYqYizyj1OpO9DfUNAsIMjf64/5gA98+H1117jpz/9KfviFF9oX4uIipKpqAAnybSJQ3C7Wi5JuV3Rx6bHEi+r3OMbF6xLqOTY62vkvn+uZ/36QGz6l7/8Je+99x6nnHKKLXvC0UIJRUkNKsCpIEKKQKHHbblpLPFKdR+Ehr27WPP7nzHypFM5aeZiXnj3c3r16tX8up2MjOApRmqooyhKYmgMOEkiTX3wNRpEAp5iIlkCqWyu4/2kkprFD2D8jfQ+9ya+OPTVxIqg/cEMC6vQhgg8PGWkiq6ipBgV4CSJNvXh4UtGJpQlMG3ikBapbbHIEQjPTDPGz95VC9j777/jLhpAUdmtuHsGju31NXLTgnXk5kjzxSPa+KIMSpRRlIxCBThJonXzSjRLIDRH2I4nHDktWPDVfEbnYd+m58TrcLnzW7xqiK+M+Y7FG9UDVpQUozHgJIlUJZbqRap4quHqd36Eb88ORITe5/2cXuf+vJX4JsKeOp9OqVCUFKMCnATla6tYWFnVqkps8pjk8mND09Ag4K3GEmFjDPvXvcjOeTez55XHA7bkuFPa9F6nVChKatEQRBJEShczkHR/Xqv99ihwc8jnb/Wa31fP7pf/l4Prl5M/cBS9zrkxqeNboVMqFCW1qAAnQbLjdMrXVjHzhY3NC2A9CtzMOH+YZdx3T52P31wyssV7GvbXsOvZO/Dt2kL3U6bSfdyliKt1Slmhx019Q0vxdruELvm51Nb5KG4qnV65qdry+Fp8oSipRQU4CZIZp1O+toppz7yLz/9VAGNPnY9pz77bqvlNkBwRykaVMHvZ5mYBduV3xZXflSMunoFn8DcsjxdaihwrKyO8wRBo8YWipAMV4CSIlC4mBEqHx81aETXtbPayzS3EN0i0zIRgdV3V7gPsq3iBLiMm4upUwJGX3h011lvocTfbYbdBetBGbbSjKOlDBTgJwtPFQj3X8LaU4SQSTy0p9FBdXU3tojvY+9EaJM9D15FnxVxoO2/EUXEfSxvtKEr60SyIJAk2uCkp9Fh2PotEtDBFjwJ3q9Jgd46w66P1HDX4OPZtWc8R595A15FntdjGSodTPbRTUZTUoAKcIuJdkJswNPKE5xyXMOP8Ydx70fGUFHoQAoJ84IM3+OhPNyOuHPpc+QDdTjiTHgVuhIBn/JtLRkYOHIfYkMjkDkVR0oeGIFKE1YJcYYGbcbNWtIqlWnmlXTvltorXjpu1gpw+x9D52G/S44wfk5PfBZ/fUJCXy9rbz2x+r1XlXHGhJ+bkDkVR2h4V4BQRaUHOnSMcONTQYox8UPQse0h4fc2C3cNXQ9+d/6aq7yRyuxXR+7xftNg2fB9WEziCkysitbi0ag6vKEr60RBEiigbVdIibFBS6KFzXm6rTIeg6FnFgINZFAc3r2Ld737CvxY9jcdbE3Hb0H3EqspLNmdZUZTUox5wCgnPHBg0fUnE7XbUenn4kpERU9j8/kZqX/sL+95+jryjvk5R2XTye5dAWAVceF5urKq8ZHKWFUVJD+oBpxErcQs+3yn3q4+/R4EbA9T880H2vf0cXUadQ5/L7iO32xHU1vlaedfhTdGjzaYbNH0JB+sbcOe0TJPQ4gpFcRb1gNNIpLiwx53DhKFFrZ4/5PMHej2MPhfP10rpMvy05tfstLaM1sTdEIgtu11CjwJ3c+mxFlcoirOoAKcRq4qy0HCBMYb9FS9Q693HwIlX02PQCQmVANtp4h4pc0JRFOdQAU4zkTzXmxasA8BfX8eXSx+hbvO/8Xz9ZGoPHuI3U8ckVAIcLvaJTmRWFKXtUAF2gOJCD1s+2kz1onto2LODwm9fRbcTJ1PSoyCpEuDQ946btUIX3RSlnaOLcA5w/Tf7suvv0/EfOsCRl95F95MupiAvN6ULYpEmHeuim6K0L9QDbkP8fj8ul4up477O+nt+w/PbOnHQ3R2AfHdqr4Xa0UxR2j8qwG1EVVUVl1xyCddddx1Tp07lxPFnsui59dC0aLanzpfy0mDtaKYo7RsNQbQBK1euZPTo0axbtw632w1ELpyI1j1NUZSOhwpwGjHGcN9993HGGWfQs2dP3n77bS6++GIg+XFGiqJkPirAaeTVV19l+vTpTJ48mbfffpvjjjuu+bVYVXKKonR8HBFgEZktIptE5D0RWSQihU7YkS72798PwIQJE3jllVdYsGABXbt2bbGNZikoiuKUB/wyMNwYcwLwH+BWh+xIOX/9618ZMGAAa9asAeC0006LODIoUve08P4OiqJ0bBzJgjDGvBTycDVwsRN2pJL6+npuvPFGHnvsMcaPH09xcXHM92iWgqJkN+0hBnw1sNTqRRG5RkQqRKSiurp9zjb77LPP+Na3vsVjjz3GLbfcwvLly+nTp4/TZimK0s5JmwcsIsuBSCr0a2PM803b/BpoAOZb7ccYMxeYC1BaWmo9s91BnnzyST744AOee+45LrzwQqfNURQlQxBjnNE0Efk+cC1wujGmzs57SktLTUVFRXoNs4nf72fbtm0MGDCAhoYGtm3bxqBBg5w2S1GUdoCIVBpjSmNt51QWxFnAL4FJdsW3PbF7927OP/98TjnlFGpra8nNzVXxVRQlbpyKAf8O6Aq8LCLrROQxh+yImwf//iLFRw/jXy++RP43vsvKTw44bZKiKBmKU1kQRztx3GS5fsYD/P7uX+Eq6E6fy++jsXgIv1q0ARHRbAZFUeKmPWRBZAR+v5/5/3iGTv2Gc9RVc+hUHCiY0P4NiqIkinZDi8GWLVtwu9307duXbufcTDd3J8TVsoJN+zcoipII6gFHYcmSJYwePZof/ehHAPQ9slcr8QXt36AoSmKoAEegsbGR22+/nfPOO4+BAwfy6KOPAtq/QVGU1KIhiDB2797N1KlTeemll/jBD37Ao48+iscT8HB1yoSiKKlEBTiM3Nxcdu7cyeOPP84Pf/jDVq9r/wZFUVKFCjCBxun/+Mc/KCsro1u3blRWVpKbqx+NoijpJetjwHV1dXz/+9/nsssuY+7cuQAqvoqitAlZrTQffvghkydPZsOGDdx555389Kc/ddokRVGyiKwV4OXLlzN58mRyc3NZunQpEydOdNokRVGyjKwNQQwYMIATTzyRNWvWqPgqiuIIWSXAX3zxBffeey/GGI455hhefvllBgwY4LRZiqJkKVkjwKtWrWL06NHceeedbNq0yWlzFEVROr4AG2N45JFHGD9+PB6Ph9WrV3Psscc6bZaiKErHF+Drr7+eG264gXPOOYeKigpGjBjhtEmKoihAFmRBTJo0iX79+nHLLbfgcnX4642iKBmEYzPhEqE9zYRTFEWxol3PhFMURVFUgBVFURxDBVhRFMUhVIAVRVEcQgVYURTFIVSAFUVRHEIFWFEUxSFUgBVFURxCBVhRFMUhMqoSTkSqga1O2xGD3kCN00akgY56XtBxz03PyzkGGGOKYm2UUQKcCYhIhZ0SxEyjo54XdNxz0/Nq/2gIQlEUxSFUgBVFURxCBTj1zHXagDTRUc8LOu656Xm1czQGrCiK4hDqASuKojiECrCiKIpDqACnGBGZLSKbROQ9EVkkIoVO25QqROS7IrJRRPwikvFpQCJylohsFpGPRGS60/akChF5UkR2icgGp21JJSLST0RWisgHTb+HNzhtU7KoAKeel4HhxpgTgP8AtzpsTyrZAFwEvO60IckiIjnAo8DZwHHAVBE5zlmrUsafgbOcNiINNAC/MMYcC4wFrsv070wFOMUYY14yxjQ0PVwN9HXSnlRijPnAGLPZaTtSxInAR8aYT4wxh4F/ABc4bFNKMMa8Dux22o5UY4z53Bizpunn/cAHQImzViWHCnB6uRpY6rQRSkRKgG0hj7eT4X/M2YSIDARGAW85a0lydPix9OlARJYDfSK89GtjzPNN2/yawC3T/La0LVnsnFsHQSI8pzmZGYCIdAEWAjcaY/Y5bU8yqAAngDHmjGivi8j3gfOA002GJVrHOrcOxHagX8jjvsAOh2xRbCIibgLiO98Y85zT9iSLhiBSjIicBfwSmGSMqXPaHsWSd4BjRGSQiOQBlwIvOGyTEgUREeCPwAfGmIecticVqACnnt8BXYGXRWSdiDzmtEGpQkQuFJHtwMnAEhFZ5rRNidK0UHo9sIzAYs7TxpiNzlqVGkTkKeBNYIiIbBeR/3LaphQxDrgSOK3pb2udiJzjtFHJoKXIiqIoDqEesKIoikOoACuKojiECrCiKIpDqAAriqI4hAqwoiiKQ2ghhuIoItILeKXpYR+gEahuenxiU5+GdouI5AI1xpgO0/VOaTs0DU1pN4jITOCAMeaBsOeFwO+q3xHDopCsAItIbkjzJiXL0BCE0i4RkaNFZENTIcsaoJ+I1Ia8fqmIPNH085Ei8pyIVIjI2yIyNsL+figiz4rIMhH5UETubXo+N8p+54nIo009aD8WkW+JyF+a+j3/MWz/D4vIGhF5ucmrR0SOaTpepYi8LiJfD9nvgyKyErgn5R+ekjGoACvtmeOAPxpjRgFVUbZ7BLjfGFMKTAGesNhuBHAxcAJwhYgU27ChuzFmAnALsBi4r8muMSIyPLgNsNoYM5pABdr/a3p+LvDfxpgxBPpC/y5kv4MJ9Aq5xYYNSgdFY8BKe+ZjY8w7NrY7g0DZbfBxDxHxGGO8Ydstb+oji4hsAvoDu2Lse3HT/+uBHcaY95ve/z4wENhEoOvdM03bzQP+3jQJZSywMMSu0L+3Z9pjSEVpW1SAlfbMwZCf/bRsIZkf8rNgb8GuPuTnRgK//9H2G/oef9j7/Xz19xO+kGKa9lljjBlpYctBi+eVLEJDEEpG0OQt7mmKq7qAC0NeXg5cF3wgIlaiF+9+7eImMKoJ4DLg38aYPcDnInJhk00uERmRwL6VDowKsJJJ/BJ4kUDa2vaQ568DxjUNQn0f+FGK9muXvcBoEVkDnArc1fT8pcC1IvIusJFAj2hFaUbT0BRFURxCPWBFURSHUAFWFEVxCBVgRVEUh1ABVhRFcQgVYEVRFIdQAVYURXEIFWBFURSH+P+GqdIHNXi7dQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x13dcb4552b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#还可以观察预测值与真值的散点图\n",
    "plt.figure(figsize=(5, 4))\n",
    "plt.scatter(y_train, y_train_pred_lr)\n",
    "plt.plot([-2.5, 2.5], [-2.5, 2.5], '--k')   #数据已经标准化，2.5倍标准差即可\n",
    "plt.axis('tight')\n",
    "plt.xlabel('True number')\n",
    "plt.ylabel('Predicted number')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "可以看出整体表现不错，相关性不太好的部分为标志线左上方的点（预测值大于真实值）和散点图最右侧的若干点（预测值小于真实值）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda\\lib\\site-packages\\sklearn\\utils\\validation.py:578: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-0.01158414,  0.30708937,  0.        , -0.03369265,  0.03610448,\n",
       "        0.00788747, -0.22574888,  0.47181581,  0.1644604 , -0.06098681,\n",
       "       -0.12121448])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 线性模型，随机梯度下降优化模型参数\n",
    "# 随机梯度下降一般不适用于数据集较小的项目\n",
    "from sklearn.linear_model import SGDRegressor\n",
    "\n",
    "# 使用默认配置初始化线\n",
    "sgdr = SGDRegressor(max_iter=1000)\n",
    "\n",
    "# 训练：参数估计\n",
    "sgdr.fit(X_train, y_train)\n",
    "\n",
    "# 预测\n",
    "sgdr_y_predict = sgdr.predict(X_test)\n",
    "\n",
    "sgdr.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of default measurement of SGDRegressor on test is -0.7711312817464651\n",
      "The value of default measurement of SGDRegressor on train is 0.7583438872459829\n"
     ]
    }
   ],
   "source": [
    "# 使用SGDRegressor模型自带的评估模块(评价准则为r2_score)，并输出评估结果\n",
    "print ('The value of default measurement of SGDRegressor on test is', sgdr.score(X_test, y_test))\n",
    "print ('The value of default measurement of SGDRegressor on train is', sgdr.score(X_train, y_train))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "性能与OLS相比差一点。因为该项目样本数较少，所以不适合使用SGD模型，sklearn建议样本数超过10万以上的大型数据集才采用SGD模型。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of RidgeCV on test is -0.6551160508536678\n",
      "The r2 score of RidgeCV on train is 0.7577250889407429\n"
     ]
    }
   ],
   "source": [
    "#岭回归／L2正则\n",
    "#class sklearn.linear_model.RidgeCV(alphas=(0.1, 1.0, 10.0), fit_intercept=True, \n",
    "#                                  normalize=False, scoring=None, cv=None, gcv_mode=None, \n",
    "#                                  store_cv_values=False)\n",
    "#class sklearn.linear_model.RidgeCV（alphas=正则参数取值范围，scoring=选择相应的评价指标，cv=交叉验证方式（留一交叉验证），gcv_mode=留一交叉验证模式）\n",
    "#导入\n",
    "from sklearn.linear_model import  RidgeCV\n",
    "\n",
    "#设置超参数（正则参数）范围\n",
    "alphas = [ 0.01, 0.1, 1, 10,100]#alphas一般在log的范围内均匀取值\n",
    "\n",
    "#生成一个RidgeCV实例\n",
    "ridge = RidgeCV(alphas=alphas, store_cv_values=True)  \n",
    "\n",
    "#模型训练\n",
    "ridge.fit(X_train, y_train)    \n",
    "\n",
    "#预测\n",
    "y_test_pred_ridge = ridge.predict(X_test)\n",
    "y_train_pred_ridge = ridge.predict(X_train)\n",
    "\n",
    "\n",
    "# 评估，使用r2_score评价模型在测试集和训练集上的性能\n",
    "print ('The r2 score of RidgeCV on test is', r2_score(y_test, y_test_pred_ridge))\n",
    "print ('The r2 score of RidgeCV on train is', r2_score(y_train, y_train_pred_ridge))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "class sklearn.linear_model.RidgeCV（alphas=正则参数取值范围，scoring=选择相应的评价指标，cv=交叉验证方式（留一交叉验证），gcv_mode=留一交叉验证模式）\n",
    "根据数据可以了解到，测试集相比于线性二乘绝对值变小，训练集比线性二乘的结果小了一些，可以推测模型加了正则函数之后出现的过拟合的情况。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl0VPeZ5vHvq40dxCKJRaxmFQYbEMSxjRc2L3F7jU3S7UzsxO1Jd9KTxJ3uuI97kpnMSafbdKeXxEnsZNLpnGRihLeQGBswxiZxgo3YLTAYsI0E2tgRm7Z3/qgrpSKXhIR0a5GezzkcV9X93ar3lkGP7u/e+15zd0RERC5VWqILEBGR1KYgERGRTlGQiIhIpyhIRESkUxQkIiLSKQoSERHpFAWJiIh0ioJEREQ6RUEiIiKdkpHoAuJh2LBhPm7cuESXISKSUjZv3nzE3XMuNq5HBMm4ceMoLi5OdBkiIinFzD5ozzhNbYmISKeEGiRmdrOZ7TGzfWb2aIzlj5jZLjPbYWbrzGxs1LIxZrbGzHYHY8YFry80sy1mts3MfmtmE8PcBhERaVtoQWJm6cATwC1AAfBJMytoMWwrUOjuM4FngMejlv0UWObu04B5QFXw+veBP3P3K4H/B/x9WNsgIiIXF+YeyTxgn7sfcPda4GngjugB7r7e3c8GTzcC+QBB4GS4+9pgXE3UOAcGBo8HAYdD3AYREbmIMA+2jwJKo56XAR9pY/xngZeCx5OBE2b2HDAeeAV41N0bgIeAVWZ2DjgFXNXVhYuISPuFuUdiMV6LeRctM7sfKASWBS9lAPOBrwBzgQnAA8GyLwO3uns+8J/At1t5z4fNrNjMiqurqy91G0RE5CLCDJIyYHTU83xiTEOZ2SLgMeB2d78Qte7WYFqsHngBmG1mOcAV7v5mMG45cHWsD3f3p9y90N0Lc3Iuehq0iIhcojCDZBMwyczGm1kW8AlgZfQAM5sFPEkkRKparDs4CA6ABcAu4DgwyMwmB68vBnaHuA0iIimp7PhZ/unld6g6fT70zwotSII9iS8Aq4n8sC9y9xIz+4aZ3R4MWwb0B1YEp/OuDNZtIDKttc7MdhKZJvth8J5/DjxrZtuBTwF/E9Y2iIikqhXFZfzg9f3UNcQ8otClQr2y3d1XAatavPa1qMeL2lh3LTAzxuvPA893YZkiIt1KY6PzzOYyrp04jFHZfUL/PF3ZLiLSzbyx/wiHTpzjvsLRFx/cBRQkIiLdTFFxGdl9M1kyPS8un6cgERHpRk6crWV1SQV3XjmKXhnpcflMBYmISDfyy22Hqa1vjNu0FihIRES6leWbSrl81EAKRg68+OAuoiAREekm3j50kl3lp1gax70RUJCIiHQbRcWlZGWkcfsVo+L6uQoSEZFu4HxdAy9sPcQtlw9nUN/MuH62gkREpBtYXVLBqfP1cT3I3kRBIiLSDRQVl5I/uA8fnTA07p+tIBERSXGlx87yxr6j3DtnNGlpse7gES4FiYhIiluxuQwz+HhhfkI+X0EiIpLCGhqdZ4pLmT8pJy4NGmNRkIiIpLA39h3h8Mnz3JegvRFQkIiIpLSi4lKy+2ayuCA+DRpjUZCIiKSo42dqWVNSGdcGjbEoSEREUtQvtx2itiG+DRpjUZCIiKQgd2d5cRkzRg2Ka4PGWBQkIiIp6O1Dp9hdfor75iZ2bwQUJCIiKamouJReGWncfsXIRJeiIBERSTXn6xp4YVvQoLFPfBs0xqIgERFJMatLKjidoAaNsShIRERSzPJNpYwe0oerEtCgMRYFiYhICik9dpbf7U9cg8ZYFCQiIilkRXFppEHjnMS1RGlJQSIikiIaGp1nNpdx3aQcRiaoQWMsChIRkRTx2+YGjclxkL2JgkREJEUUFZcyuG8miwpyE13KH1GQiIikgONnallbUsmdsxLboDEWBYmISAp4fmukQePSJGiJ0pKCREQkybk7RcWlzMwfxNThiW3QGIuCREQkye08dJJ3Kk4n3UH2JgoSEZEk19Sg8U+SoEFjLAoSEZEkdr6ugV9uO8ytM0YkRYPGWBQkIiJJ7OW3Iw0a7y1MnivZW1KQiIgkseWbShkzpC9XjU+OBo2xKEhERJLUwaNn+f2Bo9xXmJ80DRpjUZCIiCSpFZtLSTO4J4kaNMYSapCY2c1mtsfM9pnZozGWP2Jmu8xsh5mtM7OxUcvGmNkaM9sdjBkXvG5m9k0z2xss+x9hboOISCI0N2icnMOIQcnToDGW0ILEzNKBJ4BbgALgk2ZW0GLYVqDQ3WcCzwCPRy37KbDM3acB84Cq4PUHgNHA1GDZ02Ftg4hIovzm3WrKk7BBYyxh7pHMA/a5+wF3ryXyA/+O6AHuvt7dzwZPNwL5AEHgZLj72mBcTdS4vwC+4e6NwbIqRES6mRXFZQzpl8WiaXmJLuWiwgySUUBp1POy4LXWfBZ4KXg8GThhZs+Z2VYzWxbs4QBcBiw1s2Ize8nMJsV6MzN7OBhTXF1d3clNERGJn2Nnalmzq4I7rxxFVkbyH8oOs8JYpxh4zIFm9wOFwLLgpQxgPvAVYC4wgciUFkAv4Ly7FwI/BH4c6z3d/Sl3L3T3wpycnEvdBhGRuHt+6yHqGjwpGzTGEmaQlBE5ltEkHzjccpCZLQIeA2539wtR624NpsXqgReA2VHLng0ePw/MDKF2EZGEcHdWFJdyRf4gpgwfkOhy2iXMINkETDKz8WaWBXwCWBk9wMxmAU8SCZGqFusONrOmXYkFwK7g8QvBc4Drgb0h1S8iEnc7yoIGjSmyNwKRKaRQuHu9mX0BWA2kAz929xIz+wZQ7O4riUxl9QdWmBnAQXe/3d0bzOwrwDqLLNhMZBoL4B+Bn5vZl4Ea4KGwtkFEJN6KikvpnZm8DRpjCS1IANx9FbCqxWtfi3q8qI111xJj2srdTwAf68IyRUSSwrnaBlZuO8ytl49gYO/kbNAYS/KfDiAi0kO8XFLO6Qv13JsC145EU5CIiCSJ5ZtKGTu0L1dNGJLoUjpEQSIikgQ+OHqGjQeOcV/haIJjxilDQSIikgRWFJdFGjTOTu4GjbEoSEREEqypQeP1k3MYPqh3osvpMAWJiEiCbXi3mopTqdGgMRYFiYhIghVtKmVIvywWpkCDxlgUJCIiCXS05gKv7K7krlmp0aAxltSsWkSkm2hq0Jiq01qgIBERSRh3p6i4lCtGZ6dMg8ZYFCQiIgmyvewkeytrWJrCeyOgIBERSZimBo23XTEi0aV0ioJERCQBztU28Ktth7l1Rmo1aIxFQSIikgAvvR1p0JjKB9mbKEhERBJg+aZSxg3ty0fGp1aDxlgUJCIicfb+kTO8+d4x7k3BBo2xKEhEROJsxebSlG3QGIuCREQkjuobGnlmcxk3TMlNyQaNsShIRETi6DfvHqHy1AXuK+weeyOgIBERiavlm0oZ2i+LBVNTs0FjLAoSEZE46Q4NGmPpPlsiIpLknt96iPpG5765qX/tSDQFiYhIHLg7yzeVcuXobCbnpW6DxlgUJCIicbCt9ATvVtWwtJvtjYCCREQkLoqKy+iTmc5tM1O7QWMsChIRkZCdra3nV9sjDRoHpHiDxlgUJCIiIXtpZwU1F+q75bQWKEhEREK3vLiU8cP6MXfc4ESXEgoFiYhIiN47coa33jvGvYX53aJBYywKEhGREK0o7l4NGmNRkIiIhKSpQeONU3LJG9g9GjTGoiAREQnJhnerqTp9gXu7wV0Q26IgEREJyfJNpQzrn8XCabmJLiVUChIRkRAcqbnAut1V3DVrFJnp3ftHbffeOhGRBHl+S9CgsZtPa4GCRESky7k7RcWlzBqTzaRu1qAxFgWJiEgX29rUoLEH7I1AyEFiZjeb2R4z22dmj8ZY/oiZ7TKzHWa2zszGRi0bY2ZrzGx3MGZci3W/Y2Y1YdYvInIpVhSX0icznY91wwaNsYQWJGaWDjwB3AIUAJ80s4IWw7YChe4+E3gGeDxq2U+BZe4+DZgHVEW9dyGQHVbtIiKXKtKgsZyPzeyeDRpjCXOPZB6wz90PuHst8DRwR/QAd1/v7meDpxuBfIAgcDLcfW0wrqZpXBBQy4C/DbF2EZFL8uKO8m7doDGWMINkFFAa9bwseK01nwVeCh5PBk6Y2XNmttXMlgUBAvAFYKW7l3d5xSIinbSiuIwJw/pROLZ7NmiMJSPE947VncxjDjS7HygErg9eygDmA7OAg8By4AEzewm4F7jhoh9u9jDwMMCYMWM6WLqISMcdqK7hrfeP8dWbp3bbBo2xtHuPxMyuNbMHg8c5Zjb+IquUAdH7dvnA4Rjvuwh4DLjd3S9Erbs1mBarB14AZhMJlonAPjN7H+hrZvtifbi7P+Xuhe5emJOT097NFBG5ZCs2l5GeZtwzu63Jl+6nXXskZvZ1InsMU4D/BDKBnwHXtLHaJmBSEDiHgE8Af9rifWcBTwI3u3tVi3UHm1mOu1cDC4Bid38RGB61fo27T2zPNoiIhKm+oZFnN5dx45Qccrtxg8ZY2rtHchdwO3AGwN0PA21eZRPsSXwBWA3sBorcvcTMvmFmtwfDlgH9gRVmts3MVgbrNgBfAdaZ2U4i02Q/7NCWiYjE0et7e0aDxljae4yk1t3dzBzAzPq1ZyV3XwWsavHa16IeL2pj3bXAzIu8f//21CEiEramBo0LpnbvBo2xtHePpMjMngSyzezPgVfQHoKICADVpy/w6jtV3D07v9s3aIylXXsk7v7PZrYYOEXkOMnXmq7xEBHp6Z7fWhY0aOy+d0FsS3sPtvcDXnX3tWY2BZhiZpnuXhdueSIiyS3SoLGM2WOymZjb/Rs0xtLefbANQC8zG0VkWutB4CdhFSUikiq2HDzBvqqaHnUle0vtDRILWpTcDXzH3e8i0j9LRKRHK9pUSt+sdD42c2SiS0mYdgeJmX0U+DPgxeC1MK+KFxFJemcu1PPrHYf52IwR9O/Vc38ktjdIvgg8CjwXXAsyHng1vLJERJLfizvLOVPb0KOntaD9exVngUYireDvJ3KBYMy+WSIiPcWK4lIm5PRjTg9q0BhLe4Pk50SuNH+bSKCIiPRo+6tr2PT+cR69pWc1aIylvUFS7e6/CrUSEZEUsqI40qDx7h7WoDGW9gbJ183sR8A6oKlDL+7+XChViYgksfqGRp7dUsaNU3LJHdCzGjTG0t4geRCYSqTrb9PUlgMKEhHpcV7bU0316Qs99kr2ltobJFe4+4xQKxERSRHLi0sZ1r8XN/bABo2xtPf0343BfdRFRHq0qtPnefWdKu6ZPapHNmiMpb17JNcCnzaz94gcIzHA3b3NNu8iIt3N81sO0dDoPfK+I61pb5DcHGoVIiIpINKgsZQ5YwczMVe3Q2rS3jbyH4RdiIhIstty8Dj7q8/w+D2XJbqUpKIJPhGRdlre3KBxRKJLSSoKEhGRdog0aCzntpkj6NeDGzTGoiAREWmHF3eUc1YNGmNSkIiItENR0KBx9pie3aAxFgWJiMhF7KuqofiD4ywtHN3jGzTGoiAREbmIFZtLSU8z7lKDxpgUJCIibahraOTZzYdYMFUNGlujIBERacNre6o5UnOB+3Qle6sUJCIibVi+qZScAb24cUpOoktJWgoSEZFWVJ0+z/o9Vdw9exQZatDYKn0zIiKteC5o0KhprbYpSEREYnB3ijaVUjh2MJflqEFjWxQkIiIxbP7gOAeOnOE+Xcl+UQoSEZEYlm8qpV9WOh+boQaNF6MgERFpoeZCPS/uLOe2mSPVoLEdFCQiIi28uOMwZ2sbNK3VTgoSEZEWiorLuCynH7PHZCe6lJSgIBERibKv6jSbPzjO0rlq0NheChIRkSgrisvISDPumpWf6FJShoJERCRQ19DIs1vKWDA1l5wBvRJdTspQkIiIBNa/U8WRmlpdyd5BoQaJmd1sZnvMbJ+ZPRpj+SNmtsvMdpjZOjMbG7VsjJmtMbPdwZhxwes/D97zbTP7sZllhrkNItJzFBVHGjTeoAaNHRJakJhZOvAEcAtQAHzSzApaDNsKFLr7TOAZ4PGoZT8Flrn7NGAeUBW8/nNgKjAD6AM8FNY2iEjPUXXqPOv3VHPP7Hw1aOygML+tecA+dz/g7rXA08Ad0QPcfb27nw2ebgTyAYLAyXD3tcG4mqZx7r7KA8BbTeuIiHTGs80NGvUjpaPCDJJRQGnU87LgtdZ8FngpeDwZOGFmz5nZVjNbFuzhNAumtD4FvBzrzczsYTMrNrPi6urqS94IEen+3J0VxaXMGzeECWrQ2GFhBkmsE7A95kCz+4FCYFnwUgYwH/gKMBeYADzQYrXvARvc/Tex3tPdn3L3QncvzMnRfKeItK44aNB4r/ZGLkmYQVIGRJ/6kA8cbjnIzBYBjwG3u/uFqHW3BtNi9cALwOyodb4O5ACPhFS7iPQgzQ0aZ6pB46UIM0g2AZPMbLyZZQGfAFZGDzCzWcCTREKkqsW6g82saVdiAbArWOch4Cbgk+7eGGL9ItID1Fyo58Ud5fzJFSPpm6UGjZcitCAJ9iS+AKwGdgNF7l5iZt8ws9uDYcuA/sAKM9tmZiuDdRuITGutM7OdRKbJfhis8wMgD/h9sM7XwtoGEen+fr39MOfq1KCxM0KNX3dfBaxq8drXoh4vamPdtcDMGK/rVwYR6TJFxaVMzO3PrNFq0HipdLK0iPRY+6pOs+XgCZYWqkFjZyhIRKTHKmpq0Di7rSsT5GI0TdSGDXurOX62loG9MxnYJ5NBUX+yMpTBIqmsrqGR57aUsXBaLsP6q0FjZyhI2vCj377Hhr2xL2bsk5nOwD4ZfxQu0WEzsHdU8PT949d7Z6ZpN1okwdbtVoPGrqIgacO/L72SY2drOXmujpPn6jgV/DkZ48+hE+fZXX6aU+fqOH2hvs33zUpPC0InI2YIRT9vDqQgjPplpSuERLrAiuJScgf04vrJumC5sxQkbRjcL4vB/bI6vF59QyOnz9f/IYDOfzh4TjX/t54jNbXsrz7TPNZjXv8fkZFmUSGT8aEAihVKTa8N6JVBWppCSKTy1HnW76niv19/mRo0dgEFSQgy0tMuOYQaG53TF+qjgiZGCJ2v4+S5PwRV2fFzzY8bGltPITMY0Cvjj6baovd62gqlgb0z9A9Ouo1nt5TR6Ghaq4soSJJMWpo1/xDv6F9xd+dsbUPMqbfWpuUqTp7nVLD3VFvfdqOA/r0ymsNlYO+MD4XOoL6ZjB7Sl49OGErvzPQ230skUSINGsuYN34I44f1S3Q53YKCpBsxM/r1yqBfrwxGZvfp8Prn6xpiTr/FDqV6Pjh6tvm1c3UNze/TLyudG6bksmR6HjdOzWVgb917TJLHpveP896RM3z+xomJLqXbUJBIs96Z6fTOTCdvYO8Or1tb38jJc3WUHD7Jml2VrN1VyYs7y8lMN66+bBhLpuexuCCP3AEdf2+RrrR8Uyn9e2Vw64zhiS6l2zBv68huN1FYWOjFxcWJLqNHaWx0tpYeZ3VJJatLKvjg6FnMYPaYwdw0PY8lBcMZp2kFibPT5+uY98113DlrJN+6+0MdmKQFM9vs7oUXG6c9EglFWpoxZ+wQ5owdwt/dMpW9lTWsLqlgdUkF/7DqHf5h1TtMyRsQCZXpw5k+cqBOa5bQ/XpHeaRBow6ydyntkUjclR47y9pdkT2VTe8fo9FhVHYflkzP46bpw5k7bgjpOk1ZQnDX996g5nw9a758nX5xaQftkUjSGj2kL5+5djyfuXY8R2susG53FatLKvj5mwf5zzfeZ0i/LBZNy2VJwXCunTRMZ4BJl3i38jRbD57g7z82TSHSxRQkklBD+/fivrmjuW/uaM5cqOf1vdWsLqngpZ0VFBWX0TcrnRum5HDT9OHcMCWXQX10BphcmuWbSslIM+6cpQaNXU1BIkmjX68Mbp0xgltnjKC2vpGNB46yuqSCNbsqWbWzgsx046oJQ7lp+nCWFOSRewlnl0nPVFvfyPNbD7FoWp4aNIZAx0gk6UXOADvBmuBg/ftHzwIwa0w2N00fzk3Th+vCMmnTy2+X87mfbeHHDxSyYGpeostJGe09RqIgkZTi7rxbVcPqtyN7KjsPnQRgcl7/YE9lOJeP0hlg8sc+85NNlBw+yRtfXaBWPx2gg+3SLZkZk/MGMDlvAH+1cBKHTpxr3lN5Yv0+vvPqPkZl92FxQdMZYIP1g6OHqzh5ntf2VPE5NWgMjYJEUtqo7D48eM14HrxmPMfO1LJudyWrSyr5xVsH+cnv3mdw30wWTouEynydAdYjqUFj+BQk0m0M6ZfFvYWjubcwcgbYhr3VrNlVyZqSCp7ZHDkD7PrJOSyZnseCqXk6A6wHiDRoLOUj44eok0KIFCTSLfXrlcEtM0Zwy4wR1DVEnQFWUslLb1eQkWZ89LKhLAnOALuU/mKS/N567xjvHz3LXy2YlOhSujUdbJcepbHR2V52gtUlkT2VA0fOAHDl6KYzwPKYkNM/wVVKV3mkaBtrSirZ9Ngi+mRpWrOjdNZWFAWJxOLu7KuqYU3QrmVHWeQMsEm5/ZvbtcwYNUhngKWo0+frmPvNV7hrVj7funtGostJSTprS+QizIxJeQOYlDeAz984kcPNZ4BV8oPXD/DE+v2MHNS7efpr3vghOusnhfxqeznn6xpZOlcH2cOmPRKRGI6fqWXdO5EeYBv2VnOhvpHsvpksnJrHTdPzmD8pR1MlSe7OJ97gbG09q7+kBo2XSnskIp0wuF8WH5+Tz8fn5HO2tp4Ne4+wpqSCtbsqeHZLGX0y07lu8jBumj6chVPzGNRXZ4Alk72Vp9lWqgaN8aIgEbmIvlkZ3Hz5cG6+fDh1DY28eeBY0AMsMg2WkRbpAbYkuGHX8EE6AyzRlm8qJTPduEsNGuNCU1sil6ix0dlx6GTzDbsOVEfOALtidHbzXSAn5uoMsHirrW/kqm+t4yPjh/D9++ckupyUpqktkZClpRlXjs7mytHZfPXmqeyrqgmuVang8Zf38PjLe7gsp19zY8mZ+ToDLB7W7a7k2JlaXckeRwoSkS4yMbc/E3Mn8vkbJ1J+8hxrSipZs6uCJzcc4Huv7Wf4wN7NpxXPGz+ETJ0BFoqi4lKGD+zNdZNzEl1Kj6EgEQnBiEF9+PTV4/j01eM4cba2+S6QRcWl/PT3HzCoTyYLg7tAXjd5GH2z9E+xK1ScPM/re6v5yxsm6nbNcaS/vSIhy+6bxT1z8rlnTj7nahvY8G7kLpDrdlfx3JZD9MpIY/6kYSwuyGOhbrzUKU0NGu8tzE90KT2KgkQkjvpkpTcfM6lraGTTe8dYs6uStbsqeWV3FWY7mT1mMIsL8lhckMdlatfSbo2NTlFxKVdNGMLYoWrQGE86a0skCbg7u8tPs3ZX5LhKyeFTAFyW04/FBcNZXJDHrNHZpGm6plUbDxzlE09t5Nv3XcHds7VH0hV01pZICjEzCkYOpGDkQL64KHLDrleCPZUf/eYAP3h9P8P692LRtFwWF+RxzUTdW6Wlok2lDOiVwS2Xj0h0KT2OgkQkCY3K/sPB+pPn6nhtTxVrd1Xy6x3lPL2plD6ZkXurLC7IY8HUXAb3y0p0yQl16nwdq94u5+7Z+WpdkwChBomZ3Qz8O5AO/Mjd/7HF8keAh4B6oBr4jLt/ECwbA/wIGA04cKu7v29m44GngSHAFuBT7l4b5naIJNKgPpncceUo7rhyFBfqG9h44Bhrd1Xwyq4qXi6pID3NKBwbOa6ypGA4Y4b2TXTJcfer7YcjDRp17UhChHaMxMzSgb3AYqAM2AR80t13RY25EXjT3c+a2V8AN7j70mDZa8A33X2tmfUHGoNxRcBz7v60mf0A2O7u32+rFh0jke7I3dl56GTkuEpJJXsqTwMwdfiA5oP1PaUN/h3f/S3n6xp5+Uvze8T2xksyHCOZB+xz9wNBQU8DdwDNQeLu66PGbwTuD8YWABnuvjYYVxO8bsAC4E+Ddf4L+F9Am0Ei0h2ZGTPzs5mZn81fL5nCwaNnWbOrgrW7Knli/T6+8+o+hg/szaKCyPUqV00YSlZG97sI8p2KU2wvO8n/vK1AIZIgYQbJKKA06nkZ8JE2xn8WeCl4PBk4YWbPAeOBV4BHgcHACXevj3pPdWUTAcYM7ctD8yfw0PwJHDtTy6vvVEW6FW8+xM82HmRArwyunxI5rnLDlNxuc8/6ok1latCYYGEGSaxfDWLOo5nZ/UAhcH3wUgYwH5gFHASWAw8AKzvwng8DDwOMGTOmA2WLpL4hUW3wz9c18Ma+I8G1KpED9k0di5umwEZm90l0yZektr6R57eWsbggjyE9/ISDRAozSMqIHChvkg8cbjnIzBYBjwHXu/uFqHW3Rk2LvQBcBfwYyDazjGCvJOZ7Arj7U8BTEDlG0iVbJJKCemems3Ba5Kr5xkZna+mJ5imwr68s4esrS7h81EAWT4tcrzJtxICUmSJ6ZXclx8/Wca8OsidUmEGyCZgUnGV1CPgEfzi2AYCZzQKeBG5296oW6w42sxx3ryZyXKTY3d3M1gMfJ3Lm1qeBX4a4DSLdSlqaMWfsYOaMHczf3TKN/dU1rA2uV/m3dXv511f2kj+4D4um5bGkII+5Sd5csqi4lBGDenPdJDVoTKRQr2w3s1uBfyNy+u+P3f2bZvYNIqGw0sxeAWYA5cEqB9399mDdxcC/EJki2ww87O61ZjaBP5z+uxW4P2pPJiadtSVycdWnL7BudyRUfrPvCLX1jQzqk8mCqZGLIK+bnEP/Xslz6Vn5yXNc84+v8vkbJ/LXS6Ykupxuqb1nbalFioh8SNPthdfuqmTdO5WcOFtHVnoaV08MjqtMyyN3YGLvBPndV9/ln9fsZcPf3Ngjr52JBwVJFAWJyKWrb2ik+IPjzVNgB4+dBeDK0dnBRZB5TMztH9fjKo2Nzg3//Bqjsvvwi4evitvn9jQKkigKEpGu4e7sraxhbXCwfnvZSQDGDe0bnAE2nDljB4d+L5Df7z/KJ3+4kX9degV3zVKDxrAkwwWJItLNmBlThg9gyvABfGHBJCpOnmdtcFzlJ797nx/+5j2G9stqPq4yf1JOKL2viopLGdBbDRqThYJERC7Z8EG9+dS0C3WOAAAJYElEQVRVY/nUVWM5fb6O1/dWs3ZXJS+XVLBicxm9M9O4dmIOS6bnsXBqLkO74KZdJ8/VsWpnOR+fk68OyElCQSIiXWJA70xumzmS22aOpK6hkTeD5pJNF0KmGcwZO7h5Cmz8sEu7+dSvth/mQn0jS+fq2pFkoWMkIhIqd6fk8KnmO0HuLo/ctGtibn+WBFfWX5Hf/pt23f7d31Jb38hLX1SDxrDpGImIJAUz4/JRg7h81CAeWTyZ0mNneSU4rvLkhgN877X95A7oxcLgIsiPXja01Smr3eWn2FF2kq+pQWNSUZCISFyNHtKXB68Zz4PXjOfk2TrW76liza4KVm47xC/eOki/rHSumxw5rnLjlFyy+/6hh1ZRcSmZ6cadatCYVBQkIpIwg/pmcuesUdw5K3LTrt/tPxo5prKrkpfejty0a964IUHH4hxe2HqIJQXD1aAxyegYiYgkncZGZ8ehk6zdVcGakkrerappXvaTB+dyw5TcBFbXc+gYiYikrLQ048rR2Vw5Opu/uWkq7x85w9pdlRypucB8NWhMOgoSEUl644b148+vm5DoMqQVydsfWkREUoKCREREOkVBIiIinaIgERGRTlGQiIhIpyhIRESkUxQkIiLSKQoSERHplB7RIsXMqoEPLnH1YcCRLiynq6iujlFdHaO6Oqa71jXW3S/aSqBHBElnmFlxe3rNxJvq6hjV1TGqq2N6el2a2hIRkU5RkIiISKcoSC7uqUQX0ArV1TGqq2NUV8f06Lp0jERERDpFeyQiItIpCpIWzGyZmb1jZjvM7Hkzy25l3M1mtsfM9pnZo3Go614zKzGzRjNr9SwMM3vfzHaa2TYzC/22kB2oK97f1xAzW2tm7wb/HdzKuIbgu9pmZitDrKfN7TezXma2PFj+ppmNC6uWDtb1gJlVR31HD8Wprh+bWZWZvd3KcjOz/wjq3mFms5OgphvM7GTUd/W1sGsKPne0ma03s93Bv8UvxhgT7vfl7voT9QdYAmQEj/8J+KcYY9KB/cAEIAvYDhSEXNc0YArwGlDYxrj3gWFx/L4uWleCvq/HgUeDx4/G+v8YLKuJw3d00e0H/hL4QfD4E8DyJKnrAeC78fr7FPW51wGzgbdbWX4r8BJgwFXAm0lQ0w3ArxPwXY0AZgePBwB7Y/x/DPX70h5JC+6+xt3rg6cbgfwYw+YB+9z9gLvXAk8Dd4Rc12533xPmZ1yKdtYV9+8reP//Ch7/F3BnyJ/XlvZsf3S9zwALzcySoK6EcPcNwLE2htwB/NQjNgLZZjYiwTUlhLuXu/uW4PFpYDcwqsWwUL8vBUnbPkMkxVsaBZRGPS/jw//jEsWBNWa22cweTnQxgUR8X3nuXg6Rf2hAbivjeptZsZltNLOwwqY92988JvhF5iQwNKR6OlIXwD3BdMgzZjY65JraK1n/DX7UzLab2UtmNj3eHx5Mic4C3myxKNTvq0fes93MXgGGx1j0mLv/MhjzGFAP/DzWW8R4rdOnv7Wnrna4xt0Pm1kusNbM3gl+k0pkXXH/vjrwNmOC72sC8KqZ7XT3/Z2trYX2bH8o39FFtOczfwX8wt0vmNnniOw1LQi5rvZIxPd1MVuItBSpMbNbgReASfH6cDPrDzwLfMndT7VcHGOVLvu+emSQuPuitpab2aeB24CFHkwwtlAGRP9mlg8cDruudr7H4eC/VWb2PJHpi04FSRfUFffvy8wqzWyEu5cHu/BVrbxH0/d1wMxeI/LbXFcHSXu2v2lMmZllAIMIfxrlonW5+9Gopz8kctwwGYTyd6ozon94u/sqM/uemQ1z99B7cJlZJpEQ+bm7PxdjSKjfl6a2WjCzm4GvAre7+9lWhm0CJpnZeDPLInJwNLQzftrLzPqZ2YCmx0ROHIh5hkmcJeL7Wgl8Onj8aeBDe05mNtjMegWPhwHXALtCqKU92x9d78eBV1v5JSaudbWYR7+dyPx7MlgJ/LfgbKSrgJNNU5mJYmbDm45rmdk8Ij9fj7a9Vpd8rgH/F9jt7t9uZVi431e8zzBI9j/APiJziduCP01n0owEVkWNu5XI2RH7iUzxhF3XXUR+q7gAVAKrW9ZF5Oyb7cGfkmSpK0Hf11BgHfBu8N8hweuFwI+Cx1cDO4Pvayfw2RDr+dD2A98g8gsLQG9gRfD37y1gQtjfUTvr+lbwd2k7sB6YGqe6fgGUA3XB36/PAp8DPhcsN+CJoO6dtHEmYxxr+kLUd7URuDpO39W1RKapdkT93Lo1nt+XrmwXEZFO0dSWiIh0ioJEREQ6RUEiIiKdoiAREZFOUZCIiEinKEhE2mBmNZ1c/5ngqvm2xrxmbXRObu+YFuNzzOzl9o4X6QwFiUhIgl5L6e5+IN6f7e7VQLmZXRPvz5aeR0Ei0g7BFcHLzOxti9zvZWnwelrQCqPEzH5tZqvM7OPBan9G1BX1Zvb9oEFkiZn971Y+p8bM/sXMtpjZOjPLiVp8r5m9ZWZ7zWx+MH6cmf0mGL/FzK6OGv9CUINIqBQkIu1zN3AlcAWwCFgWtA+5GxgHzAAeAj4atc41wOao54+5eyEwE7jezGbG+Jx+wBZ3nw28Dnw9almGu88DvhT1ehWwOBi/FPiPqPHFwPyOb6pIx/TIpo0il+BaIl1wG4BKM3sdmBu8vsLdG4EKM1sftc4IoDrq+X1Ba/+MYFkBkbYW0RqB5cHjnwHRDfiaHm8mEl4AmcB3zexKoAGYHDW+ikirGpFQKUhE2qe1m0y1dfOpc0R6aGFm44GvAHPd/biZ/aRp2UVE9zC6EPy3gT/82/0ykR5nVxCZYTgfNb53UINIqDS1JdI+G4ClZpYeHLe4jkhzxd8SufFTmpnlEbndapPdwMTg8UDgDHAyGHdLK5+TRqT7L8CfBu/flkFAebBH9Ckit89tMpnk6P4s3Zz2SETa53kixz+2E9lL+Ft3rzCzZ4GFRH5g7yVyZ7qTwTovEgmWV9x9u5ltJdId9gDwRiufcwaYbmabg/dZepG6vgc8a2b3EunOeyZq2Y1BDSKhUvdfkU4ys/4euSveUCJ7KdcEIdOHyA/3a4JjK+15rxp3799FdW0A7nD3413xfiKt0R6JSOf92syygSzg/7h7BYC7nzOzrxO5N/bBeBYUTL99WyEi8aA9EhER6RQdbBcRkU5RkIiISKcoSEREpFMUJCIi0ikKEhER6RQFiYiIdMr/B6+2g7HUC7aWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x13dcad05c18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha is: 10.0\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'columns' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-19-06f65ed38b96>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     12\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     13\u001b[0m \u001b[1;31m# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m \u001b[0mfs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m{\u001b[0m\u001b[1;34m\"columns\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"coef_lr\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcoef_\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"coef_ridge\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mridge\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcoef_\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     15\u001b[0m \u001b[0mfs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msort_values\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mby\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'coef_lr'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mascending\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'columns' is not defined"
     ]
    }
   ],
   "source": [
    "#结果可视化\n",
    "mse_mean = np.mean(ridge.cv_values_, axis = 0)#不同alphas所对应的交叉验证的均值\n",
    "plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas),1)) \n",
    "\n",
    "#这是为了标出最佳参数的位置，不是必须\n",
    "#plt.plot(np.log10(ridge.alpha_)*np.ones(3), [0.28, 0.29, 0.30])\n",
    "\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()\n",
    "\n",
    "print ('alpha is:', ridge.alpha_)\n",
    "\n",
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef_lr\":list((lr.coef_.T)), \"coef_ridge\":list((ridge.coef_.T))})\n",
    "fs.sort_values(by=['coef_lr'],ascending=False)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "该处与前面错误相同，都为colunmns赋值时出现的问题，所以无法查看权重系数并从中进行与线性二乘的对比。\n",
    "从图中可以看出最佳alpha值为10。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of LassoCV on test is -0.8150583633359285\n",
      "The r2 score of LassoCV on train is 0.7244316704224809\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:1094: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    }
   ],
   "source": [
    "#### Lasso／L1正则\n",
    "# class sklearn.linear_model.LassoCV(eps=0.001, n_alphas=100, alphas=None, fit_intercept=True, \n",
    "#                                    normalize=False, precompute=’auto’, max_iter=1000, \n",
    "#                                    tol=0.0001, copy_X=True, cv=None, verbose=False, n_jobs=1,\n",
    "#                                    positive=False, random_state=None, selection=’cyclic’)\n",
    "from sklearn.linear_model import LassoCV\n",
    "\n",
    "#设置超参数搜索范围\n",
    "alphas = [ 0.01, 0.1, 1, 10,100]\n",
    "\n",
    "#生成一个LassoCV实例\n",
    "lasso = LassoCV(alphas=alphas)  \n",
    "#lasso = LassoCV()  \n",
    "\n",
    "#训练（内含CV）\n",
    "lasso.fit(X_train, y_train)  \n",
    "\n",
    "#测试\n",
    "y_test_pred_lasso = lasso.predict(X_test)\n",
    "y_train_pred_lasso = lasso.predict(X_train)\n",
    "\n",
    "\n",
    "# 评估，使用r2_score评价模型在测试集和训练集上的性能\n",
    "print ('The r2 score of LassoCV on test is', r2_score(y_test, y_test_pred_lasso))\n",
    "print ('The r2 score of LassoCV on train is', r2_score(y_train, y_train_pred_lasso))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LassoCV.html（LassoCV官方描述）\n",
    "alpha越大，区间越宽，可能存在的惩罚越多,该处使用的是自己设置alpha的方法。训练集测试结果相比于L2正则有一些不足"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4XVW9//H3N0Mzz0nTNk2bzjOdQgsdoEyCyCzIJCooVdSr13vFiXvFe/WqOP0QEUtFrCgCCqUigzJJCxRoU+hIS5s2adIpY5M0SYcM6/fHOS0JZDhpc7JPcj6v58nTk7NX9vnunCaf7L3WXsucc4iIiBwX4XUBIiISWhQMIiLSjoJBRETaUTCIiEg7CgYREWlHwSAiIu0ELRjM7EEzKzezzZ1sv93M1vs/NptZi5mlB6seEREJjAXrPgYzOwuoBx5yzk3tpu2lwNecc+cGpRgREQlY0M4YnHOrgOoAm18PPBKsWkREJHBRXhdgZvHARcCXA2mfmZnp8vLyglqTiMhAs27dukrnXFYgbT0PBuBS4HXnXKdnF2a2GFgMMGLECAoKCvqqNhGRAcHMdgfaNhRGJV1HN5eRnHNLnXP5zrn8rKyAAk9ERE6Sp8FgZinA2cDfvKxDRETeF7RLSWb2CLAIyDSzPcCdQDSAc26Jv9mVwPPOuYZg1SEiIj0TtGBwzl0fQJtlwLJg1SAiIj0XCn0MIiISQhQMIiLSjoJBRETaUTCIiIS4msZjLF21k7d2VfXJ64XCDW4iItKBrfvr+MPqYlas38uRplZuWzSGuaMzgv66CgYRkRDinGNNUTX3vbKTldsriI2O4MqZOXzqzDwmDU3ukxoUDCIiIaDxWDPPbTrAw2/t5u2SGjISBnH7hRP45NyRpMRH92ktCgYREY8ca27ljV1V/GPzfv6+YT/1R5sZmRHP/1w2hU/k5xI3KNKTuhQMIiJ9xDnHrsoG3tpVzeqdlax8r4JDR5uJi47ko9OG8In8XOaOSsfMPK1TwSAiEgRHm1soqWpkZ0U9W/bVsWVfHRv31FJZfxSArKQYLp42lAsmZ7NgXCax0d6cHXREwSAichJqG5tYW1xNdeMxahqPUVV/jP21RzhQd4R9NYfZV3OYVv8CmREGYwcncta4TE4flc7cUemMykzw/MygMwoGEZGTcOdTm1mxft+JzwdFRpCdEsPQ5Dhmj0zjqlnDGZUZz+jMRCYMSQqpM4LuKBhERE5Cwe6DLJqQxfcvn0pqfDSJMVEhewbQU7rzWUSkhyrrj7Ln4GHmjckgNz2epNjoARMKoGAQEemxjXtqAJg+PNXjSoJDwSAi0kPrS2uJMJiak+J1KUGhYBAR6aENpTWMz04iIWZgdtMqGEREesA5x4Y9NczIHZiXkUDBICLSIyXVjdQ0NjFdwSAiIgDrSwd2xzMEMRjM7EEzKzezzV20WWRm681si5mtDFYtIiK9ZUNpLbHREYzPTvS6lKAJ5hnDMuCizjaaWSpwH3CZc24KcE0QaxER6RXrSw8yLSeFqMiBe8ElaEfmnFsFVHfR5AZguXOuxN++PFi1iIj0hqaWVjbvqxvQl5HA2z6G8UCamb1iZuvM7FOdNTSzxWZWYGYFFRUVfViiiMj73jtwiGPNrQO64xm8DYYoYDbwMeBC4L/NbHxHDZ1zS51z+c65/KysrL6sUUTkhOMdzwN5qCp4O4neHqDSOdcANJjZKmA6sN3DmkREOrWhtIb0hEEMT4vzupSg8vKM4W/AQjOLMrN4YC6w1cN6RES6tG73QWbmpg6oCfM6ErQzBjN7BFgEZJrZHuBOIBrAObfEObfVzP4BbARagQecc50ObRUR8VL5oSPsqmzg2tNzvS4l6IIWDM656wNo81Pgp8GqQUSkt6wtOgjAnFHpHlcSfAN3IK6ISC9aU1RF/KDIATujalsKBhGRALxVVM3skWlED+Ab244b+EcoInKKahqP8V7ZIebkDfzLSKBgEBHp1trigzgXHv0LoGAQEenWmqIqBkVGDPg7no9TMIiIdGNNUTUzclOJjY70upQ+oWAQEelC/dFmNu+rC5vLSKBgEBHp0tu7D9LS6hQMIiLis6aomsgIY/bINK9L6TMKBhGRLrxVVMXUYckkxHg552jfUjCIiHSi8Vgz75TUcMaYDK9L6VMKBhGRTqwtPkhzq2P+mEyvS+lTCgYRkU6sLqwkOtLIzwuf/gVQMIiIdGr1zipmjkgjflD49C+AgkFEpEM1jcfYvK+WeWHWvwAKBhGRDr25qxrnYP7Y8OpfAAWDiEiH3thZSVx0JNOHh8f8SG0pGEREOvD6zirmjEpnUFT4/ZoMvyMWEelGed0RCsvrw7J/AYIYDGb2oJmVm9nmTrYvMrNaM1vv//husGoREemJ1TurgPDsXwAI5hisZcC9wENdtHnVOXdJEGsQEemx1TsrSYmLZtLQZK9L8UTQzhicc6uA6mDtX0QkWNYWH2TuqHQiI8zrUjzhdR/DmWa2wcyeM7MpnTUys8VmVmBmBRUVFX1Zn4iEmeaWVkqqGxmfneR1KZ7xMhjeBkY656YDvwJWdNbQObfUOZfvnMvPysrqswJFJPzsrz1CS6sjNz3O61I841kwOOfqnHP1/sfPAtFmFp49PSISMkqqGwHITY/3uBLveBYMZjbEzMz/eI6/liqv6hERgfeDYUQYB0PQRiWZ2SPAIiDTzPYAdwLRAM65JcDVwG1m1gwcBq5zzrlg1SMiEoiS6kaiIoyhKeF7KSloweCcu76b7ffiG84qIhIySqobGZ4WF7YjksD7UUkiIiGltLoxrPsXQMEgItJOaXVjWPcvgIJBROSEuiNNHGxsUjB4XYCISKgo1YgkQMEgInJCqe5hABQMIiIn6OY2n7BZ4fr1wkq++7fNTBqazORhyUwamszUYSlkJcV4XZqIhIiS6kZS4qJJiYv2uhRPhU0wxERFMDorkfWlNTy9cf+J57OTY5iWk8K0nFSm56YwfXgqaQmDPKxURLxSUn047PsXIIyCIT8vnfy8dABqDzexdX8dW/bVsXlvLZv21vLStnKO33edlxHPrJFpzBqRxul56YwbnEhEGN/sIhIuSqsbmRymazC0FTbB0FZKXDRnjM7gjNHvL9t36EgTm/bWsr60hndKali1vYLlb+8FIDU+mvyRaZw2PJWJQ5KYNDSZnNQ4hYXIANLS6thzsJELpwzxuhTPhWUwdCQpNpp5YzKZN8Y3watzjpLqRtYWH2RtUTVri6t5cWv5ifZREUZWUgyDk2NJionCDCLMMIO2ceF7zogw2rcxw3j/3w9y/hqcA4ejqcXR3NJKq4NhqbGMzkxkVGYCeZkJ5KbHERMVGcxvj8iAd6DuCE0tTpeSUDB0yswYmZHAyIwErp49HICGo81sLzvE1v2H2FvTSFndUcrqjtB4rIVW52h1vl/mxx3/pd7aCq1tfsm3tDocgPMFwPvtHf4JZ301wIkQiYowoiN9g8g27a2luqH0RLsIg2GpcYzKTGB0ZgKjsxIZO9j3MTgppt0+RaRjJVW6h+E4BUMPJMREMXNEGjNHpHldCjWNx9hV2cDuqgaKKhsprmygqLKBJ97eS/3R5hPtkmKjGJ+dxIQhSUzITmLikCQmD0smKTa8R12IfJBubnufgqGfSo0fxKwRg5j1gZByzlFx6CiFFfUUltezvewQ28vqeWbjfv58uOREu7yMeKbkpHBaTgrThqcwLSdFYSFhraS6kcgIY2hqrNeleE7BMMCYGYOTYxmcHHuivwR8gVFWd9Q/GquWLfvq2FBawzP+obtmMH5wErNGpjF7ZBpz8tLJTY/TZSgJGyXVjQxNiT1xyTacKRjChJkxJCWWISmxnDNx8InnqxuOsXFPDRtKa1lXcpCnN+7jkTW+M4shybHMGZXOmWMymDcmgxHp8QoKGbBKNKvqCQqGMJeeMIhFEwazaIIvLFpbHdvLD7G2qJo1xQd5Y1cVT23YB0BOahwLxmaycHwm88dk6kZAGVBKqxs5f1K212WEBAWDtBMRYUwckszEIcncdGYezjl2VjTwxq4qXt9RybOb9/NYQSkRBjNyUzlnwmDOmTiYKcOSdTYh/Vbt4SaqGo4xKivB61JCgoJBumRmJ4a+3nTGSJpbWtm4t5aV71Xwynvl/PyF7fz8he0MS4nl/MnZXDA5mzNGZ+g6rfQrxZUNAIzKVDBAEIPBzB4ELgHKnXNTu2h3OvAmcK1z7vFg1SO9IyoyglkjfNOFfO2C8VQcOsor75Xzwrtl/KWglIfe2E1ybBTnT8rmwqlDOHt8FrHRuvlOQluRPxhGKxiA4J4xLAPuBR7qrIGZRQJ3Af8MYh0SRFlJMVyTn8s1+bkcPtbCqzsq+OeWMl7cWsbyd/aSGBPFRyZnc+n0YSwYl6kzCQlJRZUNmMGIDHU+QxCDwTm3yszyumn2b8ATwOnBqkP6TtygSD4yZQgfmTKEppZW3txVxdMb9vPc5v0sf2cv6QmDuGz6MK6alcO0nBT1SUjIKKpsYHiappY5zrM+BjPLAa4EzkXBMOBER0awcFwWC8dl8f0rprJqewVPvrOXP68pYdnqYsZnJ/KJ/FyumJlDZqLWxBBvFVU2kJehy0jHedn5fDfwTedcS3d/OZrZYmAxwIgRI/qgNOlNg6IiOH9yNudPzqb2cBPPbNzPXwpK+cEzW/nxc9s4f1I2188dwcKxmZqxVvqcc46iygY+PivH61JChpfBkA886g+FTOBiM2t2zq34YEPn3FJgKUB+fr774HbpP1Liorlh7ghumDuC7WWH+MvaUpa/s5d/bDlATmoc18/J5bo5I3QWIX2mov4o9UebNSKpDc+CwTk36vhjM1sGPN1RKMjANT47if+6ZDK3XzSB57eU8ciaEn72/HbueamQi6cN4dPz8kJiwkIZ2Ioq/ENVsxI9riR0BHO46iPAIiDTzPYAdwLRAM65JcF6Xel/YqIiuXT6MC6dPozC8nr+9OZuHl+3hxXr95E/Mo3PLRzFBZOHEKnLTBIEGqr6YcEclXR9D9p+Jlh1SP8ydnAi37tsCl+/cAJ/WVvKg68X8YU/vc3IjHgWnzWaj88arvsipFcVVTYwKDKCYalxXpcSMjSoXEJSYkwUtywYxcrbz+G+G2eRGj+IO57czIK7/sV9rxS2W3NC5FQUVTYwMiNeZ6RtKBgkpEVGGBdPG8qKL87jz7fOZfKwZH7yj/dYcNfL/PpfCgg5dUWVDep4/gAFg/QLZsa8MZk8dMscVnxpPjNzU/npP99j4V0v88CruzjS1OJ1idIPtbQ6dlc1avK8D1AwSL8zIzeV39/sC4ipOSn84JmtnPuzV/jL2lJaWjWaWQK3r+Ywx1pa1fH8AQoG6bdm5Kbyx8/O5eHPzSUrKYZvPLGRj93zKq/uqPC6NOkndp2YVVVDVdtSMEi/N39sJiu+NJ9f3zCLhmPN3PS7NXzm92vYWVHvdWkS4or8/0fUx9CegkEGBDPjY6cN5cX/OJvvXDyRdcUHuejuVfz4uW00qINaOlFU2UBSTBSZiVqNsC0FgwwoMVGRLD5rDC9/fRGXz8hhycqdnPfzlTy3aT/Oqf9B2ttV2cCorATN9PsBCgYZkLKSYvjZNdN54rZ5pCcM4raH3+bWhwrYW3PY69IkhGioascUDDKgzR6ZxlNfns8dF0/i9cIqLvjFSpa9XkSrRi+FvSNNLeytOazptjsQcDCY2QIzu9n/OMvMRnX3NSKhICoyglvPGs0L/3EWc0al872/v8sn7n+DXeqcDmtFlQ0455uGRdoLKBjM7E7gm8C3/U9FA38KVlEiwTA8LZ7ff+Z0fnbNdLaXHeKiX77K/St36t6HMFVY7vvDQMHwYYGeMVwJXAY0ADjn9gFJwSpKJFjMjKtnD+fF/zibs8dn8aPntnH1ktUnfklI+CgsryfCNFS1I4EGwzHnG9LhAMxM30np1wYnx7L0ptn88roZFFU2cPE9r/LAq7vU9xBGCivqyU2P12y9HQg0GP5iZvcDqWZ2K/Ai8NvglSUSfGbG5TNyeOFrvrOHHzyzlRsfeIt9GrkUFnaW1zNWi/N0KKBgcM79DHgceAKYAHzXOferYBYm0leykmJYetNs7vr4NDbsqeHCu1fx9w37vC5Lgqil1bGrskH9C50ItPM5AXjZOXc7vjOFODOLDmplIn3IzLj29BE899WFjB2cyL898g7feHwDjcd01/RAVFrdyLHmVsYoGDoU6KWkVUCMmeXgu4x0M7AsWEWJeGVkRgJ/+fyZfOmcMfx13R4u+dVrvLuvzuuypJdpRFLXAg0Gc841AlcBv3LOXQlMDl5ZIt6Jjozg9gsn8vBn51J/pJkr73udx9aWaEqNAaSwQsHQlYCDwczOBG4EnvE/1+V60Wb2oJmVm9nmTrZfbmYbzWy9mRWY2YLAyxYJvnljM3n2qwvJz0vjm09s4j//qktLA0VheT2Dk2JIjtUV8Y4EGgxfBb4FLHfObfHf9fxyN1+zDLioi+0vAdOdczOAW4AHAqxFpM9kJsbw0C1z+cp543jynb1c+evVFPvn8Jf+a0d5vc4WuhBoMDQCrcD1ZrYReAo4p6svcM6tAqq72F7v3j83T8B/j4RIqImMMP7jgvEsu3kOB+qOcOm9r/HytjKvy5KT5JzzDVVVMHQq0GB4GHgQXx/DpcAl/n9PiZldaWbb8F2euuVU9ycSTGePz+Lpf1tAblo8n/1DAfe8tEP9Dv1QWd1R6o82Kxi6EGgwVDjn/u6cK3LO7T7+caov7px70jk3EbgC+H5n7cxssb8foqCiQss2indy0+N54rZ5XDEjh1+8sJ0v//kdDh9r8bos6YETI5J0c1unuuxAbuNOM3sAX7/A0eNPOueW90YRzrlVZjbGzDKdc5UdbF8KLAXIz8/Xn2jiqbhBkfziE9OZOCSJH/9jG8VVDfz2U/kMS43zujQJQGH5IUAjkroS6BnDzcAMfJ3Jl/L+5aSTZmZjzb9skpnNAgYBVaeyT5G+YmZ8/uwxPPjp0ympauSye19nQ2mN12VJAAor6kmKjSIrKcbrUkJWoGcM051z03qyYzN7BFgEZJrZHuBOfNN145xbAnwc+JSZNQGHgWudLthKP3POxMEs/+I8bl62lmuXvsHd187goqlDvS5LulDo73jWcp6dC/SM4U0z69ENbc65651zQ51z0c654c653znnlvhDAefcXc65Kc65Gc65M51zr/W4epEQMC47iRVfms+kocl84U9v85tXdqpTOoQVljeof6EbgQbDAmC9mb3nvyltk3/Yqojgu9/hkVvP4JLThnLXP7bxnSc309zS6nVZ8gG1jU1U1h9V/0I3Ar2U1NWNaiICxEZHcs91MxmRHs99r+xkX81hfn3jLBJjAv0xk2ArrPB1PI/LVjB0JdBpt3d39BHs4kT6m4gI4xsXTeTHV03jtcJKrlnyBmV1R7wuS/x2lB0fqqoFKLsS6KUkEemB6+aM4MHPnE5JVQNX3aelQ0NFYXk9sdER5KRpaHFXFAwiQXL2+CweXXwmR5tbuHrJatbtPuh1SWFvR3k9ozMTiYzQiKSuKBhEgmja8BSW3zaf1Lhobvjtm7y0VXMseamwvF79CwFQMIgE2YgM3zQaE4YksfiP6/hrQanXJYWlhqPN7K05rKGqAVAwiPSBjMQY/nzrGZw5OoPbH9/I/St3el1S2NnpX5xHZwzdUzCI9JHEmCh+95l8LjltKD96bhs/enarboTrQ+8v56kRSd3RAGuRPhQTFckvr5tJanw096/axcHGY/zwymlERepvtGDbUV5PVIQxMiPe61JCnoJBpI9FRhjfv3wq6Qkx3PPSDmoPN/HL62YSGx3pdWkDWmF5PaMyE4hWCHdL3yERD5j5VoW789LJ/HNLGbcsW0v9Ua0nHUyFWrUtYAoGEQ/dPH8U/+/a6bxVVM0Nv32Tqvqj3X+R9NiRphZ2VzUwTsEQEAWDiMeunDmcpTfN5r0Dh7jm/jfYW3PY65IGnOKqBlodjFEwBETBIBICzpuUzR8/O5eKQ0e5+jer2VF2yOuSBpTjcySN04ikgCgYRELEnFHp/OXzZ9Lc6rh6yRuaQqMXFZbXYwajsxK8LqVfUDCIhJBJQ5NZfts80uKjufEBTaHRWwrL6xmRHq+RXwFSMIiEmNz0eB6/bR7jBidx60MFPLqmxOuS+r0d5Yc0FUYPKBhEQlBmYgyPLj6DheOy+NbyTdz94nbdJX2SmltaKapsYKymwgiYgkEkRCXERPHAp/O5evZw7n5xB7c/vpFjzVoutKd2VzfS1OJ0xtADQQsGM3vQzMrNbHMn22/0rx+90cxWm9n0YNUi0l9FR0bw06tP46vnjePxdXu46XdvUdN4zOuy+pWt++sAX/+NBCaYZwzL6Hqt6CLgbOfcacD3gaVBrEWk3zIzvnbBeO6+dgbvlNRw5X2r2VWhFeECtXV/HVERprueeyBoweCcWwVUd7F9tXPu+Hi8N4HhwapFZCC4YmYOD986l9rDTVz+69f517Zyr0vqF97dV8eYrESNSOqBUOlj+CzwXGcbzWyxmRWYWUFFRUUfliUSWk7PS+dvX5pPblo8t/xhLfe+vEOd0t3Yuv8Qk4bqxrae8DwYzOwcfMHwzc7aOOeWOufynXP5WVlZfVecSAjKTfetCHfZ9GH87PntfO4PBZpjqRMHG45xoO6I+hd6yNNgMLPTgAeAy51zVV7WItKfxA2K5O5rZ3DnpZN5dUclH/3lq7y2o9LrskLO8Y7nycMUDD3hWTCY2QhgOXCTc267V3WI9Fdmxs3zR7HiS/NJio3ik797i//5+xYaNH33Ce9qRNJJCeZw1UeAN4AJZrbHzD5rZl8wsy/4m3wXyADuM7P1ZlYQrFpEBrLJw5J5+t8WctMZI/n968Vc8IuVvPCuptIAXzBkJcWQmRjjdSn9StBWcHPOXd/N9s8BnwvW64uEk7hBkXz/iqlcMXMY31m+mVsfKuDciYP5+kcmhPVlFF/Hc/ge/8nyvPNZRHrP7JHpPP2VBXz7oxMpKK7m4nte5ct/fpvC8vCbxvtYcyuF5YeYrGDoMa35LDLAREdG8Pmzx3DdnBH8dtUuHny9iKc37mfemAw+ecZILpicHRbrHheW19PU4jRU9SQoGEQGqJS4aL5+4QRunp/HYwWlPPxmCV98+G2SY6OYMyqDM0anc8HkbEZmDMw1Ck6MSNIZQ48pGEQGuIzEGL64aCyfP2sMK7eX8/yWMt7cVcWLW8v4+fPbeeK2eQOyH2Lr/jpioiIYlTkwgy+YFAwiYSIywjh3YjbnTswGoLiygeuWvsniPxbw1JcXkJ4wyOMKe9fWA3VMGJJEVBhcNutt+o6JhKm8zATuv2k25YeO8uU/v01zy8CZ0ts5x7v76pg0ZOCdCfUFBYNIGJuem8oPr5zG6p1V/N+zW70up9eU1R3lYGOTOp5PkoJBJMxdPXs4n5mXx+9fL2Z14cCYVmPLvloAJg9L8biS/knBICJ866MTycuI544VmznS1OJ1OadsQ2kNEQZTc3Qp6WQoGESE2OhIfnDFNIoqG7jvX4Vel3PK1u+pZXx2EvGDNL7mZCgYRASABeMyuXJmDr9ZubNf3yntnGPjnhqmD0/1upR+S8EgIifc8bFJJMRE8Z3lm2lt7Z8LAJVUN1LT2MT0XAXDyVIwiMgJmYkxfOejk1hTXM1f15V6Xc5JWV9aA8D0XHU8nywFg4i0c03+cOaMSueHz26jsh+uDLehtJbY6AjGZ2uo6slSMIhIO2bGD6+cSuOxZn7w9Ltel9NjG/fUMGVYSlhMFBgs+s6JyIeMHZzEbYvGsmL9Pl7dUeF1OQFramll875adTyfIgWDiHToi4vGMDozgTue7D/3NmwvO8SRplb1L5wiBYOIdCg2OpIfXDmVkupGfvXyDq/LCcjGPb47nnXGcGoUDCLSqXljMvn4rOHcv3IX28tC/96GDaU1pMRFMzIj3utS+rWgBYOZPWhm5Wa2uZPtE83sDTM7amZfD1YdInJq7vjYJJJio/jO8k0hf2/D+tIapuemYmZel9KvBfOMYRlwURfbq4GvAD8LYg0icorSEwbxnYsnUbD7II8VhO69DY3HmtlRXs+M4epfOFVBCwbn3Cp8v/w7217unFsLNAWrBhHpHVfPHs7cUen86NmtVBwKzXsbtuyro6XVcZr6F06Z+hhEpFtmxg+vmsbhphZ+FKLrNqwt9v0dOnOEguFU9YtgMLPFZlZgZgUVFf1nTLXIQDImK5HPnzWG5e/s5c1dVV6X8yFri6oZk5VARmKM16X0e/0iGJxzS51z+c65/KysLK/LEQlbXz53LLnpcfz3is0caw6dpUBbWh0Fuw8yZ1S616UMCP0iGEQkNMRGR/K9S6ewo7ye371W5HU5J7x34BCHjjRzep6CoTcEbRULM3sEWARkmtke4E4gGsA5t8TMhgAFQDLQamb/Dkx2ztUFqyYROXXnTcrmI5OzueelHVxy2lBy072/Z+B4/4KCoXcEc1TS9c65oc65aOfccOfc75xzS5xzS/zbD/ifT3bOpfofKxRE+oE7L5tCZITxtcfW09zi/SWlNcXVDE2JZXhanNelDAi6lCQiPZaTGscPrphKwe6D3OvxUqDOOdYWVXN6XrpubOslCgYROSlXzMzhqpk53PPSDgqKO71lKehKqhspP3SU09Xx3GsUDCJy0v7n8ikMT4vnq4+up/awN/eqrinyhdIc9S/0GgWDiJy0pNhofnndDMrqjvCVR96hxYO5lNYWV5MSF824wYl9/toDlYJBRE7JzBFp/O/lU1m5vcKTu6LXFh/k9Lw0IiLUv9BbFAwicspumDuCz8zL44HXinhsbUmfvW75oSMUVTZomGovUzCISK/4r49NYsHYTP5rxeYT1/2DraD4IIA6nnuZgkFEekVUZAS/vmEWuWnxfOFP6yitbgz6a75WWEliTBTTcjTVdm9SMIhIr0mJj+aBT+fT3NLKZ/+wlkNHgjtS6bUdlZwxOp3oSP0q6036bopIrxqdlch9N85mZ0UDX310fdBGKpVUNVJS3cjCcZpYs7cpGESk1y0Yl8n3Lp3My9vK+ck/twXlNV4trDhG6Ge3AAALFklEQVTxWtK7gjaJnoiEt5vOzGPbgUPcv3IX03JSuOS0Yb26/9d2VDIsJZbRmQm9ul/RGYOIBNGdl05h1ohUvvH4Rt47cKjX9tvS6li9s4oF4zI1P1IQKBhEJGgGRUXwm0/OJiEmisV/LKC2sXc6ozftraX2cBML1L8QFAoGEQmq7ORYfnPjLPYePMx//nUDzp16Z/Sr2339C/PHZJzyvuTDFAwiEnT5eel8++JJvLi1jD+sLj7l/b1aWMnUnGSt7xwkCgYR6RO3zM/jvImD+eGz29i8t/ak99NwtJl3Sg6yYKwuIwWLgkFE+oSZ8dNrppOWEM1XHnmHhqPNJ7Wft4qqaGpxLNQw1aBRMIhIn0lPGMQvr5tJcVUD33tqy0nt45X3KoiNjmD2yLRerk6OC1owmNmDZlZuZps72W5mdo+ZFZrZRjObFaxaRCR0nDE6g9sWjeGv6/bw0tayHn2tc46XtpazYGwWsdGRQapQgnnGsAy4qIvtHwXG+T8WA78JYi0iEkK+et54Jg5J4lvLN3Gw4VjAX7e9rJ69NYc5f9LgIFYnQQsG59wqoKu5dy8HHnI+bwKpZjY0WPWISOgYFBXBLz4xg5rGY3y3B5eUXvSfYZw7UcEQTF72MeQApW0+3+N/TkTCwORhyXzl3HH8fcM+ntm4P6CveXlbOacNT2FwcmyQqwtvXgZDR/exd3jni5ktNrMCMyuoqKgIclki0lduWzSG6cNT+K8VmyivO9Jl26r6o7xdclBnC33Ay2DYA+S2+Xw4sK+jhs65pc65fOdcflaWxi6LDBRRkRH84toZHG5q4fbHN3Z5V/Qr71XgHJw3MbsPKwxPXgbDU8Cn/KOTzgBqnXOBnU+KyIAxJiuROy6exMrtFfzxzd2dtntpWxnZyTFMzUnuw+rCU9Cm3TazR4BFQKaZ7QHuBKIBnHNLgGeBi4FCoBG4OVi1iEho++QZI3lpWzn/98xW5o3JYOzgpHbbjzW3smp7JZdOH6rZVPtA0ILBOXd9N9sd8KVgvb6I9B9mxk8+fhoX3r2Krz66nuVfnEdM1Pv3Kawtrqb+aLMuI/UR3fksIiFhcHIsP7l6Olv21fHj595f9e1IUwu/fHEHcdGRzB+raTD6goJBRELGBZOz+cy8PH7/ejEvvFtGS6vj3x9dz5riau66+jTiBulu576gpT1FJKR8++KJrC2u5vbHN7BofBb/2HKA714ymcum9+7SoNI5nTGISEiJiYrk3htm0dTcyor1+/jC2WO4ZcEor8sKKzpjEJGQMyozgSU3zWbT3lpuO3uM1+WEHQWDiISkheOyWKg1nT2hS0kiItKOgkFERNpRMIiISDsKBhERaUfBICIi7SgYRESkHQWDiIi0o2AQEZF2rKsVk0KRmVUAna/m0bsygco+eq1gGijHATqWUDVQjmWgHAd8+FhGOucCumOw3wVDXzKzAudcvtd1nKqBchygYwlVA+VYBspxwKkdiy4liYhIOwoGERFpR8HQtaVeF9BLBspxgI4lVA2UYxkoxwGncCzqYxARkXZ0xiAiIu0oGNows++b2UYzW29mz5tZh2sJmtmnzWyH/+PTfV1nd8zsp2a2zX8sT5pZaiftis1sk/94C/q6zkD04FguMrP3zKzQzL7V13UGwsyuMbMtZtZqZp2OFukn70ugxxLS74uZpZvZC/6f5RfMLK2Tdi3+92O9mT3V13V2pbvvsZnFmNlj/u1vmVletzt1zunD/wEkt3n8FWBJB23SgV3+f9P8j9O8rv0DNX4EiPI/vgu4q5N2xUCm1/We6rEAkcBOYDQwCNgATPa69g7qnARMAF4B8rto1x/el26PpT+8L8BPgG/5H3+ri5+Veq9rPdnvMfDF47/LgOuAx7rbr84Y2nDO1bX5NAHoqAPmQuAF51y1c+4g8AJwUV/UFyjn3PPOuWb/p28Cw72s51QEeCxzgELn3C7n3DHgUeDyvqoxUM65rc6597yuozcEeCz94X25HPiD//EfgCs8rOVkBPI9bnuMjwPnmZl1tVMFwweY2f+ZWSlwI/DdDprkAKVtPt/jfy5U3QI818k2BzxvZuvMbHEf1nSyOjuW/vaedKe/vS+d6Q/vS7Zzbj+A/9/BnbSLNbMCM3vTzEIpPAL5Hp9o4/8jqxbI6GqnYbfms5m9CAzpYNMdzrm/OefuAO4ws28DXwbu/OAuOvjaPh/a1d1x+NvcATQDD3eym/nOuX1mNhh4wcy2OedWBafizvXCsYTEewKBHUsA+s370t0uOngupH5WerCbEf73ZDTwspltcs7t7J0KT0kg3+Mevw9hFwzOufMDbPpn4Bk+HAx7gEVtPh+O7zprn+ruOPyd4pcA5zn/xcUO9rHP/2+5mT2J77S0z38B9cKx7AFy23w+HNjXexUGrgf/v7raR794XwIQEu9LV8dhZmVmNtQ5t9/MhgLlnezj+Huyy8xeAWbiu7bvtUC+x8fb7DGzKCAFqO5qp7qU1IaZjWvz6WXAtg6a/RP4iJml+UcwfMT/XMgws4uAbwKXOecaO2mTYGZJxx/jO47NfVdlYAI5FmAtMM7MRpnZIHwdbCE1ciRQ/eV9CVB/eF+eAo6PLPw08KEzIf/Peoz/cSYwH3i3zyrsWiDf47bHeDXwcmd/LJ7gda96KH0AT+D7IdwI/B3I8T+fDzzQpt0tQKH/42av6+7gOArxXVNc7/84PiJhGPCs//FofCMYNgBb8F0e8Lz2kzkW/+cXA9vx/RUXqsdyJb6/3o4CZcA/+/H70u2x9If3Bd+19peAHf5/0/3Pn/iZB+YBm/zvySbgs17X/YFj+ND3GPhffH9MAcQCf/X/LK0BRne3T935LCIi7ehSkoiItKNgEBGRdhQMIiLSjoJBRETaUTCIiEg7CgYJG2ZWf4pf/7j/zteu2rzS1Wyjgbb5QPssM/tHoO1FTpWCQSQAZjYFiHTO7err13bOVQD7zWx+X7+2hCcFg4Qd8/mpmW32r3twrf/5CDO7z7/OwNNm9qyZXe3/shtpc1esmf3GP6naFjP7n05ep97Mfm5mb5vZS2aW1WbzNWa2xsy2m9lCf/s8M3vV3/5tM5vXpv0Kfw0iQadgkHB0FTADmA6cD/zUP0/OVUAeMA34HHBmm6+ZD6xr8/kdzrl84DTgbDM7rYPXSQDeds7NAlbSft6tKOfcHODf2zxfDlzgb38tcE+b9gXAwp4fqkjPhd0keiLAAuAR51wLUGZmK4HT/c//1TnXChwws3+1+ZqhQEWbzz/hnxI7yr9tMr6pVNpqBR7zP/4TsLzNtuOP1+ELI4Bo4F4zmwG0AOPbtC/HN92ESNApGCQcdbZISVeLlxzGN+cMZjYK+DpwunPuoJktO76tG23nnznq/7eF938Ov4Zv3qHp+M7mj7RpH+uvQSTodClJwtEq4Fozi/Rf9z8L3+RirwEf9/c1ZNN+evWtwFj/42SgAaj1t/toJ68TgW82S4Ab/PvvSgqw33/GchO+ZRuPG0//nWVV+hmdMUg4ehJf/8EGfH/Ff8M5d8DMngDOw/cLeDvwFr7VrsC3Nsci4EXn3AYzewff7Ke7gNc7eZ0GYIqZrfPv59pu6roPeMLMrgH+5f/6487x1yASdJpdVaQNM0t0ztWbWQa+s4j5/tCIw/fLer6/byKQfdU75xJ7qa5VwOXOt864SFDpjEGkvafNLBUYBHzfOXcAwDl32MzuxLd+bklfFuS/3PULhYL0FZ0xiIhIO+p8FhGRdhQMIiLSjoJBRETaUTCIiEg7CgYREWlHwSAiIu38f7/IpNdt4VX4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x13dcb046a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha is: 0.16703196638946255\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'columns' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-48-5b33142c6968>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      9\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m \u001b[1;31m# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 11\u001b[1;33m \u001b[0mfs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m{\u001b[0m\u001b[1;34m\"columns\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"coef_lr\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcoef_\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"coef_ridge\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mridge\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcoef_\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"coef_lasso\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlasso\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcoef_\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     12\u001b[0m \u001b[0mfs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msort_values\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mby\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'coef_lr'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mascending\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'columns' is not defined"
     ]
    }
   ],
   "source": [
    "mses = np.mean(lasso.mse_path_, axis = 1)\n",
    "plt.plot(np.log10(lasso.alphas_), mses) \n",
    "#plt.plot(np.log10(lasso.alphas_)*np.ones(3), [0.3, 0.4, 1.0])\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()    \n",
    "            \n",
    "print ('alpha is:', lasso.alpha_)\n",
    "\n",
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef_lr\":list((lr.coef_.T)), \"coef_ridge\":list((ridge.coef_.T)), \"coef_lasso\":list((lasso.coef_.T))})\n",
    "fs.sort_values(by=['coef_lr'],ascending=False)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "alpha的值较大，虽然因为columns的原因看不到权重系数，但是可以根据alpha值估计，可能会出现权重系数为0的情况。"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "整个作业的完成过程比较困难，问题不断，但是通过查询交流等方式也对所学知识有了拓展和补充。因为最近家中有事，学习不太系统，作业完成也比较仓促，没有过多的加入新的东西，只是对所学知识进行了回顾和重现，希望助教谅解。\n",
    "\n",
    "作业中存在如下问题：\n",
    "1.column赋值报错问题如何解决\n",
    "2.当r2_score出现负值时应该如何判断模型的优劣\n",
    "3.如果先切分数据如何与已经定义好的X和y进行整合\n",
    "谢谢助教，辛苦了！"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
